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ABSTRACT 


Anode  sheaths  impact  the  operation  of  many  practical  plasma  devices.  This 
complex  region  is  explored  in  detail  for  collisional.  isothermal  (identical  specie 
temperatures),  low-temperature  plasmas,  where  sheath  dimensions  are  in  the  micron 
range.  The  selected  approach  involves  postulation  of  a  specific  electric  field 
distribution  with  two  shape  factors.  Previous  research  regarding  planar  anodes  is 
verified  and  expanded  upon  using  greater  parameter  ranges,  'z',  a  dimensionless 
quantity  specifying  plasma  composition  and  condition,  groups  diverse  plasmas  into 
■families'  exhibiting  similar  sheath  characteristics,  'rj'.  a  nondimensional  ratio  of 
electrical  energy  to  thermal  energy  in  the  sheath,  allows  temperature  effects  to  be 
studied.  The  investigation  focuses  on  three  disparate  plasma  families  that  span  a  z 
range  of  1.1729  to  2.1493,  at  q  values  defined  by  plasma  temperatures  of  bOOO'  K, 
3000"  K.  and  300'  K.  Results  indicate  that  at  lower  temperatures,  charge  production 
in  the  outer  sheatli  is  generic  to  the  electric  field  distribution,  and  that  the  sheaths 
themselves  are  nearly  unaffected  by  substantial  changes  in  temperature  (i.e.,  q). 
Conversely,  sheath  density  and  extent  are  shown  to  vary  significantly  for  differing  z 
values.  Newly-derived  equations  governing  cylindrical  anodes  generate  sheaths  that 
are  virtually  identical  to  corresponding  planar  cases.  It  is  shown  that  only  those 
anodes  whose  radii  are  comparable  to  the  plasma's  'characteristic  radius'  (y)  must  be 
treated  with  the  cylindrical  formulation;  non-vacuous  plasmas  would  require  micron- 
width  anodes  to  be  thus  affected.  Finally,  an  analytical  approach  yields  solutions  that 
confinn  the  numerical  results,  and  offers  an  algebraic  approximation  for  high-q 
plasmas. 
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I.  INTRODUCTION 


Various  types  of  plasma  devices  have  been  in  operation  for  over  30  years. 
Marshall's  coaxial  plasma  gun  [Ref.  1]  was  successfully  accelerating  volumes  of 
hydrogen  plasma  as  far  back  as  the  late  1 950's.  In  the  intervening  years,  the  number 
of  applications  for  plasma  devices  has  steadily  increased  to  include 
magnetohydrodynamic  (MHD)  power  generation,  laser  pumping,  strategic  defense, 
and  electromagnetic  propulsion  for  interplanetary  spacecraft.  However,  due  to  their 
somewhat  complex  nature,  comparatively  little  is  understood  about  the  sometimes- 
destructive  sheath  regions  surrounding  the  electrodes  in  every  plasma  device.  This 
work  attempts  to  further  previous  research  efforts  concerning  description  of  the 
anode  sheath. 

As  stated  by  Biblarz  [Ref.  2],  completely  satisfactory  anode  sheath  solutions  do 
not  exist  for  several  plasma  conditions:  one  such  case  involves  steady,  collisional, 
low-temperature,  isothermal  discharges.  He  then  goes  on  to  derive  an  involved, 
nonlinear  differential  equation  that  describes  the  entire  plasma  region  affected  by  a 
planar  anode,  from  the  surface  to  the  undisturbed  plasma.  A  presumed  (but 
'shapeable')  function  describing  the  electric  field  in  that  region,  selected  after  much 
deliberation  [Ref.  2  and  Ref.  3],  ultimately  allows  charge  production  rates  and 
electron/ion  populations  to  be  plotted  as  a  function  of  distance  from  the  anode. 
Nondimensional  parameters  make  the  solution  profiles  applicable  to  several  families 
of  plasmas. 

In  this  work,  numerical  techniques  are  employed  to  verify  the  one  previously 
solved  case,  and  to  explore  planar  anode  solutions  to  several  other  plasma  conditions. 
This  can  be  readily  accomplished  due  to  the  generality  of  the  formulation  and  the  ease 
with  which  the-profiles  can  be  produced.  Sheaths  and  ambipolar  regions  are  properly 


generated  for  all  cases.  The  effect  of  nondimensional  parameter  variations  on  the 
sheath  is  extensively  investigated,  and  some  general  conclusions  are  hypothesized.  In 
addition,  an  analytical  technique  is  presented  that  supports  the  numerical  results  and 
allows  for  an  accurate  algebraic  solution  to  low-temperature  conditions.  Finally, 
similar  derivations  produce  equations  that  treat  the  cylindrical  anode  sheath  problem: 
the  equations  and  their  profiles  are  then  analyzed,  and  a  comparison  is  made  to  planar 
anode  findings. 


II.  THE  CARTESIAN  PROBLEM 


This  section  discusses  and  validates  one  previous  anode  sheath  research  effort 
concerning  steady,  low-temperature  collisional  plasmas.  In  addition,  several  new 
parameter  cases  are  examined  and  analyzed.  The  geometry  of  the  Cartesian  problem 
is  illustrated  in  Figure  1. 


Figure  1  Cartesian  Anode  Geometry 


A.  DERIVATION  REVIEW 

Previous  work  by  Biblarz  [Ref.  2]  treats  the  planar  anode  sheath  problem  in  a 
one-dimensional  Cartesian  fashion,  and  is  reviewed  here  as  it  forms  the  starting  point 
for  this  work.  Electric  field  intensity  (E),  electron  and  ion  population  densities  (  n , 
and  nj  ),  and  all  other  relevant  quantities  are  considered  to  vary  only  with  linear 
distance  (y)  from  the  flat  anode  surface.  The  applicable  relations  are  Gauss'  equation 
and  the  two  species  continuity  equations,  as  shown  on  the  following  page. 
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Note  that  e  is  the  electron  charge  constant.  £„  is  the  permittivity  constant,  the 

subscripts  e  and  i  indicate  electron  and  ion  particles,  the  j’s  are  the  species 
contributions  to  the  total  current  density  (i.e.,  J=j,.+j,),  the  p's  are  the  particle 

mobilities,  and  the  D's  are  the  diffusion  coefficients.  Although  this  set  does  not  fully 
constrain  the  problem  (as  an  equation  describing  plasma  reactivity,  h,.,  is  absent),  any 

selection  of  a  specific  form  for  E(y)  implicitly  fixes  that  variable.  It  is  of  interest  to 
explore  the  problem  with  'guesstimates’  of  E(y)  in  order  to  facilitate  solutions. 

Reference  2  then  introduces  two  new  parameters,  designated  K*  and  K  ,  which 
are  defined  by  Equations  (2a),  (2b),  and  (2c).  These  K-terms  are  somewhat  artificial 
parameters,  although  they  are  related  to  the  total  current.  More  meaningful  is  the 
derivative  of  K",  which  is  directly  proportional  to  the  net  production  rate  of 
charges,  h^  (i.e.,(K*)  ~t\JD^  ). 
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Manipulation  of  the  previous  six  equations  and  the  Einstein  relation  (which 
relates  the  mobilities  to  the  diffusion  coefficients)  yields  the  isothermal  differential 
equation  below.  Note  that  k  is  the  Boltzmann  constant,  and  that  all  primes  denote 
derivatives  with  respect  to  the  variable  y. 


.  2J 
+  K  = - + 

filLl 

EE- 

rkT;,£,YE"'  E'E"Y 

U. 

eD 

r 

1  i  E  E^  ;j 

Advanced  knowledge  of  E(y)  considerably  facilitates  the  solution  of  Equation  (3). 
Biblarz'  justification  for  the  selected  form  of  E(y)  is  detailed  elsewhere  [Ref  2  and 
Ref.  3].  This  function  meets  several  critical  conditions  inherent  to  the  stated 
problem,  including  the  required  distribution  of  E  and  the  behavior  of  both  K-terms  at 
the  boundaries.  The  function  is  reproduced  as  Equation  (4)  below. 


E(y)  =  E. 


(4) 


The  parameters  ’a’  and  ’A'  are  shape  factors,  where  'a'  is  of  the  order  of  the 
sheath  length,  and  'A'  relates  to  the  physical  parameters  A4)  and  E,  (as  discussed  on 

pages  7  and  8). 

Substitution  of  the  chosen  E(y)  into  Equation  (3).  followed  by  an  order-of- 
magnitude  analysis,  the  nondimensionalization  of  the  parameters  in  accordance  with 
Figure  2.  and  further  simplification  based  on  field  and  current  properties  at  the 
electrode  and  in  the  undisturbed  plasma,  all  combine  to  yield  the  governing 
differential  equation  and  boundary  condition  of  Figure  3. 
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Figure  2  Nondimensionalized  Parameters 
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Figure  3  Simplified  Nondimensional  Differential  Equation 
and  Boundary  Condition  Governing  the  1-D  Cartesian  Problem 


The  introduction  of  the  dimensionless  parameters  n  and  z  further  generalizes  the 
problem,  in  that  one  solution  to  the  above  equation  can  be  app.icable  to  a  large 
number  of  specific  dimensional  cases.  Their  additional  significance  is  discussed  in  a 
subsequent  section. 

As  will  be  seen,  both  numerical  and  analytical  solutions  to  the  equations  of  Figure 
3  are  possible,  yielding  profiles  of  K*  and  (K’)  thai  are  a  function  of  distance  from 
the  planar  anode. 

Further  manipulation  of  Equations  (1)  and  (2)  derives  relations  that  Produce 
electron  and  ion  population  profiles  [Ref.  2].  A  more  precise  fonn  of  these  equations 
is  given  in  Figure  4.  Note  that  specific  solutions  for  K*  are  required  to  generate  die 
h -profiles.  Several  distinct  cases  in  the  next  section  yield  population  curves  that 
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clearly  illustrate  tlie  propriety  of  the  sheatli  and  ambipolar  regions  in  the  plasma  near 
the  anode. 


Reference  2  also  presents  two  more  useful  relations  from  analysis  involving 
physical  observations  (Equations  (5)  and  (6)).  Recent  analysis  has  proven  that  the 
infinite  series  in  Equation  (6b)  converges  for  ah  values  of?.  These  are  combined  and 
manipulated  to  yield  tlie  two  important  equations  of  Figure  5.  which  allou 
detennination  of  the  parameters  z  and  q.  The  procedure  is  as  follows: 

•  the  type  of  gas  defines  the  anode  potential  drop  A(|)  (which  is  essentiall) 

equivalent  to  the  ioniz.ation  potential  of  tlie  gas) 

•  tlie  particular  case  or  application  specifies  E.  and  n 

•  z  is  then  found  by  iteration  of  tlie  implicit  equation  at  the  top  of  Figure  5 

•  the  choice  of  X  locks  in  q  using  the  bottom  equation  in  Figure  5 


f  2eE  ) 


V  en  , 


b  exp(2z 


A(j)  =  H  -Za  f(z) 


where  f(z)  =  z2]7 


z 


‘-,(2m  -  l)ni: 


(5) 

(6a) 

(6b) 
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Figure  5  Equations  That  Yield  /  and  r| 


Specific  values  of  z  and  q  allow  computation  of  the  electric  field  parameters  A 
and  a  (Equations  (5)  and  (6)).  which  in  turn  yield  the  tailored  form  of  the  E 

distribution  (Equation  4),  and  ultimately  permit  the  generation  of  K*.  [K  )  .  and  h^ 
profiles  in  the  plasma  regions  close  to  the  anode  (equations  of  Figure  3  and  Figure  4). 

B.  NUMERICALLY  SOLVED  CASES 

1.  Procedure 

Numerical  solution  of  the  governing  differential  equation  (Figure  3)  can  be 
achieved  with  a  fourth-order  Rungc-Kutta  FORTRAN  program.  First  the  dependent 
variable  is  redefined  and  the  equations  are  rewritten  in  accordance  with  Figure  6. 


^[w]  =  2q  -q 
d> 


vE„y 


w-  q 


(E/EJ 

(y  +  1)" 


w(0)  =  1 


where  w(y)  =  — - — 

^  E(y)/E., 

Figure  6  Modified  1-D  Cartesian  Differential  Equation 
and  Boundary  Condition  for  Numerical  Solution 
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Application  of  the  Runge-Kutta  scheme  computes  the  value  of  w  at  each  \ . 
The  data  for  K*  is  then  recovered  with  the  division  of  each  w  datapoint  by  the 

corresponding  value  for  E(y  )/E„.  (K  )  data  are  extracted  with  the  following  steps 
for  each  value  of  y ; 

•  compute  w*  using  the  defining  differential  equation  (Figure  6) 

•  perform  the  operation  below  (derived  from  the  Quotient  Rule) 


dv 


1 

(  eY 

f  E'l 

' 

/E  ^ 

w' 

—  +  K 

IeJ 

1  eJ 

(7) 


Finally,  the  h^,  profiles  are  computed  as  a  function  of  y  using  the  available 

information  and  the  equations  of  Figure  4.  A  single  computer  routine  can  be  made  to 
perfonn  all  of  the  required  operations  (see  Appendix  A);  the  resulting  data  are  then 
plotted. 

Note  again  that  this  technique  is  substantially  simpler  than  alternate  methods 
tliat  do  not  presume  a  specific  form  of  E(y).  Ensuing  discussions  address  the 
advantages  and  disadvantages  of  this  approach. 

2.  Case  I:  the  Nitrogen  Problem 

The  first  specific  case  involves  verification  of  an  earlier  example  [Ref  2].  A 
planar  anode  in  contact  with  nitrogen  plasma  was  analyzed  under  the  following 
conditions  (nitrogen  is  a  common  discharge  gas  with  well-known  properties): 

•  nitrogen's  anode  potential  drop  (Acj),):  15.51  V  (singly  ionized) 

•  the  electric  field  strength  in  the  undisturbed  plasma  (E. ):  12.(XX).  V/m 

•  the  particle  density  (n, ):  ICf' m  ‘ 

•  the  temperature:  bOOG^K 
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Tlie  expression  used  to  approximate  f(z)  (Equation  6a)  was  tamcated  after  13 
terms,  producing  an  error  in  the  calculation  of  z  which  is  on  the  order  of  10  . 

Using  the  given  data,  the  following  parameters  result: 

•  z  =  1.75626 

•  r|  =  99.12285 

•  a  =  1.95421  10'^  m  and  A  =  1.17793  x  10  "  m 
.  E„  =  262,265.  V/m 

The  computed  value  for  E,,  (obtained  using  Equation  4)  relates  specifically  to 
all  plasmas  that  are  defined  by  the  above  values  of  z  and  q;  this  includes  nitrogen 
plasma  at  the  stated  conditions.  The  normalized  electric  field  E(y  )  for  this  case  is 
shown  in  Figure  7.  Note  that  the  field  decreases  monotonically  and  abruptly  to  a 
constant  value  in  the  undisturbed  plasma,  as  is  required  of  the  model. 

I 

The  corresponding  K  ,  (K  )  .  and  h^.  ^  profiles  are  presented  in  Figure  8  and 
Figure  9.  Note  that  those  Figures  also  represent  all  other  plasma  cases  whose  gas 
composition,  density,  temperature,  and  electric  field  intensity  are  defined  by 
z=  1.75626  and  n=99. 12285.  Figure  9  clearly  illustrates  the  distinct  sheath  and 
ambipolar  regions  present  for  this  general  case.*  The  K-term  profiles  of  Figure  8 

likewise  exhibit  expected  characteristics,  with  K*  rising  roughly  monotonically  from 

/ 

+  1  to  +2.  while  (K  )  decreases  to  near-zero  at  the  outer  edge  of  the  sheath.  The 

leftmost  downturn  of  (K* )  ,  while  possibly  correct,  is  quantitatively  suspect  since  the 
model  of  ionization  by  electron  impact  breaks  down  at  the  fringe  of  the  collisionless 
region  (i.e.,  for  y/a  s  -10  ). 


*  The  sheath  is  that  small  region,  extending  away  from  the  anode,  where  the 
electron  and  ion  populations  are  not  equivalent;  in  the  ambipolar  region,  the  rauo 
of  oppositely-charg^  particles  is  1:1,  but  the  populations  of  both  are  less  than  those 
that  exist  in  the  'undisturbed  plasma’ 
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(lu/A)  (»VX)3 


Fifjure  7  Normalized  Electric  Field  (z=  1.75626,  ri=^91229) 


Figure  8  K-term  profiles  (7=1.75626,  r|=99.1229) 
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Figure  9  n, ,  profiles  (z=l. 75626,  r|=99.1229) 
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As  the  parameter  'a'  specific  to  Case  1  is  known,  it  can  be  seen  that  the 
nitrogen  plasma  sheath  for  the  particularly  defined  discharge  extends  approximately 
0.1  mm  out  from  the  planar  anode.  Similarly,  the  region  of  undisturbed  plasma  for 
such  a  device  begins  approximately  6  mm  from  the  anode  surface.  Both  of  these 
results  are  in  accordance  with  expectations,  the  generalities  of  which  were  obtained 
after  much  deliberation  [Ref.  4]. 

These  numerically-obtained  Case  I  results  confirm  the  work  of  Reference  2. 
Postulation  of  a  suitable  ionization/recombination  mechanism  is  left  to  that  discourse. 

The  remaining  objectives  of  this  work  are  as  follows: 

•  investigation  of  the  effects  which  are  produced  as  the  defining  parameters  are 

varied  over  a  wide  yet  practical  range 

•  presentation  of  an  analytical  solution  to  the  general  Cartesian  anode  sheath 

problem 

•  derivation  and  examination  of  the  related  Cylindrical  anode  sheatli  problem 

3.  Cases  II  and  III:  Varying  Temperature  in  the  Nitrogen  Problem 

Comprehensive  investigation  of  the  nitrogen  plasma  problem  requires 
consideration  of  temperature's  effect  on  the  anode  sheath.  Thus,  while  all  other  Case 
I  conditions  arc  held  constant,  two  somewhat  more  practical  values  for  the  isothermal 
temperature  are  considered: 

•  Case  II  X,  “t  3000"K 

•  Case  Ill  sets  X,  at  300"K  (this  is  representative  of  discharges  in  laser  pumping) 

Temperature-only  variations  impact  r|  exclusively;  ?  and  the  normalized 
electric  field  E(y  )/E„  (including  the  inherent  factors  A  and  a)  are  unaffected.  Since 
q  is  by  definition  inversely  proportional  to  X  (Figure  2).  the  decreased  temperatures 
of  Cases  II  and  III  signify  larger  values  of  q-  As  a  result,  the  first-order  derivative 
term  m  the  governing  differential  equation  (Figure  3)  becomes  less  influential  at  the 
lower  temperatures.  Indeed,  for  the  conditions  of  Case  III.  the  contribution  of  the 
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derivative  term  is  negligible  (note,  however,  that  this  circumstance  does  not  nullify 
the  validity  of  the  formulation:  the  equation  does  not  become  degenerate). 

Application  of  the  lower  equation  in  Figure  5  yields  q  values  of  198.246  and 
1982.46  for  Case  II  and  Case  III,  respectively.  The  K-tenn  profiles  for  these  general 
cases  are  presented  in  Figures  10  and  1 1.  Both  cases  exhibit  the  same  general  K-term 
tendencies  that  have  been  recognized  previously,  with  minor  variations  that  illustrate 
the  decreasing  importance  of  the  net-ionization  tenn. 

The  effect  of  temperature  variation  (or,  equivalently,  q  variation)  on  a 
plasma  in  the  anode  region  is  explored  using  Figures  12  and  13.  which  directly 

compare  the  K  and  (K  )  profiles  of  the  three  cases.  One  significant  observation, 
contnbuted  by  Biblarz  subsequent  to  Reference  2,  hypothesizes  that  the  asymptotic 
behavior  of  the  K-tenns  with  decreasing  temperature  reflects  the  approach  to 
ionization  that  is  genenc  to  the  electric  field  distribution;  the  overall  results  are  thus 
independent  of  the  details  of  ionization  and  recombination.  As  an  important  practical 
example,  measurements  of  charge  density  would  not  reflect  relevant  details  of  the 
ioniz.ation/recombination  mechanisms. 

Conversely,  the  h.  profiles  for  this  general  z=1.7.‘>626  case  are  nearly 
unaffected  by  the  changes  in  temperature.  Figure  14  illustrates  charged  particle 
curves  that  exhibit  only  a  barely  perceptible  shift  of  position  over  the  entire  range  of 
temperatures  tested.  This  unexpected  result  has  substantial  implications;  the  most 
significant  of  these  is  the  acceptability  of  the  isothermal -particle  assumption  (i.e.,  the 
assertion  that  X  =T  =X.)  Though  such  a  premise  may  be  unrealistic,  its  validity 

appears  to  be  inconsequential  to  the  population  profiles.  Electric-field  effects 
dominate  any  temperature  effects  in  the  cases  examined.  In  fact,  the  physical 
significance  of  q  (stated  in  Reference  2  as  a  ratio  of  electrical  to  thermal  energy) 
predicts  such  dominance  with  its  large  magnitudes. 


I.S 


Fijjure  10  K-terni  profiles  (/.=  !. 75626,  r|=198.246) 
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Kplus  Profil 


Figure  12  K*  profiles  (/,=  1. 75626) 
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Figure  13  (K  )  profiles  (/.=  1.75626) 
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Figure  14  n^.  profiles  (z=l. 75626) 
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C.  EXPLORING  THE  ENVELOPE  FOR  z  AND  n 

The  foregoing  results  include  a  preliminary  investigation  into  the  effect  of  p- 
variation  on  the  planar  anode  sheath.  This  upcoming  section  more  fully  explores  how 
changes  in  the  parameters  z  and  q  impact  the  sheath  region. 

1.  Practical  Extremes  of  z 

The  parameter  z  is  a  dimensionless  representation  of  three  basic  variables 
(A4)^  ,  E.  ,  and  nj  that  specify  the  condition  of  a  plasma;  several  different  gases,  at 

differing  densities  and  field  strengths,  can  be  related  by  their  equivalent  z  values. 
Various  diverse  plasmas  belonging  to  the  same  'z-family'.  at  the  same  temperature, 
would  theoretically  exhibit  identical  sheath  characteristics.  Using  the  technique 
discussed  previously,  anode  sheaths  for  plasmas  defined  by  two  limiting  values  of  z 
are  examined. 

Singly-ionized  elements  possess  anode  potential  drops  (A(t)J  that  range  from 

24.46  Volts  (helium)  down  to  3.87  Volts  (cesium),  with  the  low  values  being 
preferred.  Common  densities  (n.)  for  these  type  collisional  plasmas  vary 

approximately  from  10  '  to  lO  "  particles  per  cubic  meter.  Finally,  reasonable  field 
strengths  in  the  undisturbed  plasma  (E. )  are  somew'hat  arbitrarily  set  from  5000  to 

30.IX)0  Volts  per  meter.  Application  of  the  upper  equation  in  Figure  5  yields  the 
following  'practical'  extremes  for  z; 

1.1729  s  z  s  2.1493  (8) 

This  relatively  small  numerical  range  represents  an  extremely  diverse 
grouping  of  plasmas;  z  is  a  decidedly  sensitive  parameter.  Small  z  values  represent 
those  plasmas  w  ith  small  potential  drops  and  low  densities  that  are  subjected  to  very' 
high  electric  fields  (e.g..  low'  density  cesium  at  30.0(X).  V/m).  Plasmas  with  large  z 
values  are  characterized  by  gases  w'ith  large  potential  drops,  at  conditions  of  high 
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densin  and  lou  field  strength  (e.g.,  high  pressure  helium  at  5000.  V/m).  By  way  of 
comparison,  the  nitrogen  plasma  considered  earlier  in  this  work  (median  potential 
drop  and  density,  small  field  strength)  has  a  slightly -larger-than-median  z  value  of 
1.7563. 

Values  for  p  are  required  to  completely  constrain  the  problem,  which  entails 
the  choosing  of  one  or  more  temperatures.  To  facilitate  direct  comparisons,  the  same 
arbitrary  quantities  for  X  that  have  been  considered  for  the  nitrogen  case  are 

designated  as  the  standard  values  (i.e.,  6000" K,  3000" K.  and  300'  K). 

Consequently,  numerical  sheath  solutions  are  next  generated  and  examined 

for  these  forementioned  limiting  conditions: 

•  'Small  z'  plasmas  (z=  1.1729)  at  the  three  standard  temperatures  (q=  16.4135. 

32.8269,  and  328.269  for  X  at  6000"K,  3000"K,  and  300" K) 

•  'Large  z'  plasmas  (z=  2.1493)  at  the  standard  temperatures  (r)=  257.563. 

5 1 5. 1 25,  and  5 1 5 1 .25  as  X  decreases  to  300"K ) 

2.  Sheath  Solutions  for  'Small  z'  Plasmas 

K-term  profiles  for  the  'small  z'  plasmas  (at  the  three  standard  temperatures) 
are  depicted  in  Figures  15.  16.  and  17.  Qualitatively,  each  of  these  curves  displays 
the  same  previously  noted  characteristics  inherent  in  the  nitrogen  plasma  curves; 

monotonically  increasing  and  well-behaved  K’.  and  a  rate-of-charge-production  term 

/ 

(K")  that  decreases  and  vanishes  at  the  outer  edge  of  the  sheath.  As  before,  the 

leftmost  downturn  of  (K  )  is  quantitatively  suspect  at  the  fringe  of  the  collisionless 
region. 

These  reassuring  results  extend  to  the  composite  curves  of  Figures  18  and  19. 
which  depict  K-term  profile  variations  as  a  function  of  temperature  (i.e.,  q)  changes. 
In  particular,  the  charge  production  rate  curves  of  Figure  19  again  illustrate  the 
striking  implication  that  charge  production  in  the  outer  sheath  is  independent  of  the 


Figure  15  K-term  profiles  {z=1.1729,  q  =  I6.4135) 
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Figure  16  K-term  profiles  (z=1.1729,  ti=32.8269) 
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K-term  Profiles,  Low  z,  T=300°J 


TT-n 

10 


Figure  17  K-terni  profiles  (z=1.1729,  ri=328.269) 


Figure  18  K'  profiles  (z=1.1729) 
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(Kplus)  Profiles,  Low  z,  Various  Temp. 


Figure  19  (K‘)  profiles  (z=1.1729) 
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temperature.  Moreover,  the  inner  sheath  once  again  sees  an  increase  in  the  net 
production  of  charges  as  the  temperature  decreases. 

Electron  and  ion  population  profiles  for  the  'small  z'  plasmas  are  given  in 
Figure  20.  As  with  the  nitrogen  case  (and  other  plasmas  of  that  z-family).  the  sheath 
and  ambipolar  regions  are  clearly  visible  and  quite  appropriate,  while  the  population 
curves  themselves  again  appear  to  be  nearly  unaffected  by  large  variations  in 
temperature. 

Such  qualitative  observations  for  the  'small  z'  plasmas  indicate  that  previous 
impressions  concerning  'median  z'  plasmas  (e.g.,  nitrogen)  are  not  limited  to  that  one 
category,  but  are  possibly  characteristic  of  aU  collisional  low-temperature  plasmas. 
The  sheath  profiles  generated  for  the  'large  z'  cases  of  the  next  section  further 
confirm  this  hypothesis. 

Quantitative  deviations  in  the  sheath  due  to  the  change  of  z  are  addressed  in  a 
later  section. 

3.  Sheath  Solutions  for  'Large  z'  Plasmas 

The  high  density,  small  field  'large  z'  plasmas  continue  the  trends  established 
by  the  other  collisional  plasmas,  albeit  with  some  minor  differences.  One  such 
difference  is  illustrated  by  the  K-term  profiles  of  Figures  21.  22.  and  23.  The  curs  es 
display  the  same  general  traits  of  the  other  plasmas,  with  the  exception  that  the  K* 
curve  no  longer  climbs  monotonically.  This  small  anomaly  becomes  more 
pronounced  at  higher  temperatures  (Figure  21).  yet  does  not  appear  to  significantly 
impact  the  charge  production  rate  curves,  which  still  behave  as  expected. 

Similarly,  there  are  no  surprises  in  the  composite  curs'es  of  Figures  24  and 

I 

25.  In  fact,  the  (K  )  curves  of  Figure  25  demonstrate  even  more  strikingly  the 
generic,  electric-field-dependent  nature  of  ionizabon  and  recombination  in  the  sheath. 
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n-Profiles,  Low  z,  Varying  Temp 


Figure  20  j  profiles  (z=1.1729) 


Figure  21  K-lerm  profiles  (/=2.1493,  n=257.563) 
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Figure  23  K-term  profiles  (/=2.1493,  q =5 15 1.2 5) 


Figure  24  K*  profiles  (z=2.1493) 


Fijiure  25  (K  )  profiles  (7=2.1493) 
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The  n-profiles  of  tliis  class  of  plasmas  are  depicted  in  Figure  26.  While  there 
are  some  yet-to-be-addressed  differences  from  previous  n-profiles.  the  sheath  and 
ambipolar  regions  are  again  evident,  as  are  the  profiles’  theorized  independence  from 
temperature  changes. 

4.  The  Effect  of  z  on  the  Planar  Anode  Sheath 

Results  to  this  point  have  revealed  some  important  and  perhaps  non-intuitive 
concepts  concerning  the  properties  of  planar  anode  shec-Jis.  concepts  which  appear 
applicable  to  ah  low-temperature,  collisional  plasmas.  However,  sheath  variations 
attributable  solely  to  changes  in  the  parameter  z  have  yet  to  be  discussed. 

Comparison  of  the  n-profiles  for  each  z-family  (Figures  14,  20.  and  26) 
illustrates  a  subtle  point:  although  all  n-profiles  have  been  shown  to  be  nearly 
independent  of  temperature,  they  do  become  slightly  more  sensitive  to  temperature 
variations  as  the  value  of  z  decreases.  This  is  possibly  attributable,  in  part,  to  the 
lower  densities  of  'small  z'  plasmas;  the  temperature-dependent  kinetic  energy 
changes  of  the  particles  may  produce  a  more  observable  effect  on  the  total 
iomzation/recombination  process  in  a  plasma  not  already  saturated  with  density- 
driven  collisions.  Perhaps  more  contributory  is  the  previously-discussed  dominance 
of  electrical  energy  over  thermal  effects.  The  'small  z'  plasmas  exhibit  smaller 
magnitudes  of  E„  than  other  plasma  families,  decreasing  such  dominance.  The 
involvement  of  A4)^  in  this  phenomenon  remains  unclear.  In  any  event,  such  n- 
profile  temperature  shifts  remain  nearly  negligible  for  all  plasmas  considered  to  this 
point,  which  further  validates  the  innocuous  nature  of  the  constant  temperature 
assumption. 

Figure  27  depicts  the  K*  curves  of  three  diverse  z-families  at  the  same 
temperature  (6(X)0‘’K).  As  noted  previously,  all  are  well-behaved  and  do  not  differ  to 
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Figure  26  fi,,  profiles  (z=2.1493) 
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any  significant  degree.  Contrastingly,  the  charge  production  rate  curves  of  these 
three  plasma  families  (Figure  28)  are  somewhat  disparate  and  very  insightful.  'Large 
z'  plasmas  show  significant  charge  production  throughout  the  entire  sheath,  in 
contrast  to  the  localized  and  decreased  production  rates  depicted  for  the  'small  z' 
plasmas.  This  may  indicate  the  fact  that  electric-field  effects  are  more  localized  for 
smaller  z. 

The  n-profiles  for  the  three  z-families  are  depicted  in  Figure  29.  Unlike 
similar  profiles  for  varying  temperature  (i.e.,  q).  these  curves  are  decidedly  affected 
by  changes  in  z.  Most  prominently,  the  magnitude  of  the  charged  particles  changes 
dramatically  as  z  varies,  especially  in  the  sheath.  This  is  a  graphic  indication  of  z- 
induced  changes  in  both  E„  and  current  flow  for  diverse  plasmas.  Also  noteworthy  is 

the  fact  that  the  sheath  itself  extends  further  from  the  anode  surface  as  the  value  of  z 
decreases.  In  particular,  the  sheath  of  the  'small  z'  plasmas  stretches  over  ten  times 
the  distance  occupied  by  the  'large  z'  sheath. 

5.  Some  Final  Thoughts  on  the  Influence  of  q 

The  effect  of  q  variations  on  the  anode  sheath  has  already  been  indirectly 
addressed  via  the  extensive  consideration  of  temperature's  influence  on  the  problem. 
The  bottom  equation  of  Figure  5  inversely  relates  these  two  parameters  for  fixed  z 
and  A({), .  Physical  interpretation  of  q  as  a  ratio  of  electrical  and  thermal  energy  has 

also  been  reviewed. 

The  purely  numerical  importance  of  q  is  apparent  with  a  glance  at  the 
governing  differential  equation  (Figure  3);  its  reciprocal  controls  the  influence  of  the 
charge  production  rate  term.  Large  values  of  q  reduce  this  term  to  near  zero,  the 
consequences  of  which  can  be  seen  in  every  low-temperature  case  (i.e.,  300‘’K).  To 
fully  and  completely  explore  this  parameter's  effect  on  the  sheath,  solutions  are 
presented  for  the  original  nitrogen  plasma  at  artificially  small  values  of  q  (in 
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n-Profiles  as  a  Function  of  z,  T=600CfK 


Figure  29  n,.  Profiles  at  6000" K  for  Different  z-Families 


particular,  n  is  set  equal  to  25.,  10..  and  1.0).  Note,  however,  that  tlie  validity  of 

these  curves  is  questionable  at  best,  since  such  small  q  values  equate  to  some 

extraordinarily  high  temperatures  for  the  nitrogen  (up  to  595.000  K!).  the  existence 

of  which  contradicts  the  original  assumption  of  low -temperature  plasmas. 

Nevertheless,  the  following  curves  may  perhaps  reveal  some  valid  trends. 

Despite  the  unusual  and  probably  erroneous  K*  plots  for  the  small-q 
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conditions  (Figure  30).  the  (K  )  curves  of  Figure  31  hearteningly  retain  some  of  the 
familiar  traits  concerning  charge  production  rate  in  the  sheath.  The  corresponding  n- 
profiles  are  presented  in  Figure  32;  remarkably,  only  the  anomalous  'ri=r  profiles 
exhibit  tendencies  somewhat  contrary  to  those  already  obser\  ed. 

D.  ANALYTICAL  SOLUTIONS 

In  order  to  confirm  the  validity  of  the  numerically-obtained  solutions,  an 
analytical  solution  to  the  problem's  governing  relation  (Figure  6)  is  presented. 

1.  Procedure  for  'The  Outer  Expansion  .Method' 

Reference  5  offers  an  approach  which  is  ideally  suited  to  the  form  of  the 
differential  equation  in  Figure  6.  The  term  [l/q)  is  justly  defined  as  a  'small' 

parameter,  and  thj  e(,;’ation  is  already  in  dimensionless  fonn.  A  series  solution  for  w 
is  presumed,  approximate  yet  valid  for  'non-small'  values  of  y.  Equation  (9)  depicts 
the  power  senes  expansion  for  w : 


r  1 3 

f\\ 

w  =  1  w  j 

=  W„  + 

- 

|w  +  j 

w.  + 

InJ 

.nj 

.n, 

Substitution  of  Equation  (9)  into  the  governing  equation  allows  tenns  with 
like  powers  of  (l/n)  to  be  equated,  which  in  turn  allows  each  w^  to  be  solved  for 
analytically.  Because  the  solution  lor  w„  is  purely  algebraic  in  this  fonnulation.  the 
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Figure  30  K*  profiles  (/=  1.75626,  Small  n*') 
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n-Profiles,  z  for  Nitrogen,  Low  eta’s 


Figure  32  profiles  (z=l. 75626,  wSmall 


first-order  derivatives  in  each  of  the  subsequent  \v.  expressions  can  be  evaluated 
exactly  (i.e.,  symbolically).  The  total  solution  w  is  thus  an  infinite  summation  of 
exact  expressions.  Prevalent  magnitudes  of  r|  allow  highly  precise  approximations  of 
w  with  the  series  truncated  to  only  three  terms. 

Utilization  of  this  technique  produces  the  following  approximate  analytical 
'outside'  solution  (w„ )  to  the  governing  differential  equation  of  Figure  6: 
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Figure  33  Three-Term  Approximate  Analytical  Solution  (Wy) 


The  structure  of  the  equations  themselves  offers  insight.  Note  that  each  w’^  is 
a  function  only  of  z  and  y ;  they  are  completely  independent  of  r|  (i.e.,  temperature). 

The  contribution  of  q  to  the  total  solution  manifests  itself  solely  in  the  series- 
expanded  expression  for  w' .  Note  also  that  the  equation  for  w„  is  identical  to  the 

original  governing  differential  equation  for  the  case  when  q  is  allowed  to  increase  to 
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infinity.  As  a  result,  the  partial  solution  \v„  is  itself  a  valid  approximation  for  the 

true  solution  (w)  at  large  values  of  q*  This  fact  is  also  corroborated  by  the  series- 
expanded  expression  for  w/.  all  terms  beside  w„  vanish  for  large  q. 

As  before,  plots  of  K*  can  be  recovered  from  the  data  for  w  (or,  in  this  case, 
w_).  Using  parameter  values  from  the  first  nitrogen  plasma  case  (z=  1.7563  and 

q=99.1229),  K"  profiles  from  both  the  'Outer  Solution'  and  the  Fourth-Order  Runge- 
Kutta  scheme  are  compared  (Figure  34).  Although  the  ’(Xiter  Solution'  does  indeed 
appear  to  diverge  for  extremely  small  values  of  y,  the  two  curves  are  nearly 
identical.  Such  correlation  offers  comforting  evidence  for  the  validity  of  the 
numerical  procedures  and  results  presented  in  this  work. 

Figure  35  illustrates  the  soundness  of  using  just  w,,  to  recover  data  for  large- 

q  cases.  The  solid  curve  is  from  the  Runge-Kutta  solution  for  the  300" K  nitrogen 
case  (z=1.7563  and  q=1982.5);  the  dotted  profile  depicts  data  generated  using  only 
the  partial  solution  w„. 
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Figure  34  Comparison  of  Numerical  and  Analytical  Solutions 
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Figure  35  One  Term  (w„)  Algebraic  Approximation  for  High  r|*^ 
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III.  THE  CYLINDRICAL  PROBLEM 


The  following  sections  present  detailed  derivations  and  analysis  for  the 
cylindrical  anode  sheath  problem.  As  with  the  planar  cases,  the  focus  is  limited  to 
steady,  low-temperature  collisional  plasmas.  Exploration  of  the  cylindrical  sheath  is 
desirable,  as  several  functioning  plasma  devices  contain  or  employ  a  cylindrical 
architecture.  In  addition,  these  results  can  be  compared  and  contrasted  to  the 
previously-presented  Cartesian  results.  The  geometry  of  the  one-dimensional 
cylindrical  problem  is  illustrated  in  Figure  36. 


A.  DERIVATIONS 

Attainment  of  a  governing  differential  equation  for  the  cylindrical  case  begins 
with  the  re-derivation  of  the  applicable  relations  (Gauss'  equation  and  the  two  species 
continuity  equations)  in  cylindrical  form.  Reference  6  presents  the  following 
differential  form  of  Gauss'  law : 


V.E  =  i 

c . 


(10) 
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V  is  tlie  del  operator,  E  is  the  vector  form  of  the  electric  field,  p  is  the  charge 
density  (which  can  be  represented  here  in  accordance  with  Equation  1 1 ).  and  £,  is  the 

permittivity  of  free  space.  The  divergence  (V  E)  of  the  electric  field  in  cylindrical 
coordinates  is  given  in  Equation  (12)  below; 

p  =  e(n  -nj  (11) 


V  E  = 


c 


\rjCT 


1  VE  fE 


I  \  1  L  I  L 


ae 


cz 


(12) 


For  a  one-dimensional  electric  field  that  varies  only  as  a  function  of  radius  (r). 
combination  of  Equations  (10).  (11),  and  (12)  yields  the  following  one-dimensional 
cylindrical  fonn  of  Gauss’  equation 


dE 

dr 


(13) 


Comparison  of  Exjuations  (13)  and  (la)  accentuates  the  appearance  of  a  new  term 
which  is  the  result  of  the  cylindrical  derivation. 

The  species  continuity  equations  are  derived  from  conservation  equations  [Ref.  7] 
that  can  be  manipulated  into  the  form  given  below; 

-j,  = -ep,n,E-eD^Vn^.  (14a) 

j=ep,nE-eDVn  (14b) 

For  a  one-dimensional  electric  field,  the  previous  relations  are  easily  rewritten  as 
show'n  in  Equations  (15a)  and  (15b)  on  the  following  page;  these  are  the  1-D 
cylindrical  species  continuity  equations. 
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j  =e|j  n  E  +  eD  — ^ 
dr 


(15a) 


dn 

j  =e|j  n  E-eD 

dr 


(15b) 


Note  that  these  species  equations  are  nearly  identical  to  those  for  tlie  planar  case 
(Equations  lb  and  Ic);  only  the  independent  variable  has  been  renamed. 

The  procedure  from  this  point  is  the  same  as  that  of  Reference  2: 

•  Combine  these  forms  of  Gauss'  equation  and  the  species  continuity  equations 

with  the  Einstein  relation 

•  Incorporate  the  K-terms  as  defined  by  Equation  (2) 

•  Assume  isothermal  conditions 

•  Derive  a  single  nonlinear  differential  equation  in  K',  E.  and  their  derivatives 
The  desired  intermediate  differential  equation  (the  cylindrical  counterpart  to 

Equation  3)  is  presented  below.  Note  that  all  primes  here  denote  derivatives  with 
respect  to  tlie  variable  r. 


1  .  2J  i 

^  £  Y 
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'  E"  ^  E'  ' 
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.rE;  Vr-'Ey 
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f  ^ 


(16) 


In  comparison  to  Equation  (3),  several  new  terms  exist  that  are  due  solely  to  the 
cylindrical  nature  of  the  derivation.  However,  the  validity  of  the  above  expression  is 
demonstrated  by  letting  r  increase  towards  infinity,  which  causes  the  equation  to 
approach  the  planar  case.  In  that  limit.  Equation  (16)  reduces  exactly  to  the 
corresponding  Cartesian  expression  of  Equation  (3). 

The  next  step  calls  for  an  'educated  guess'  of  the  form  for  E(r).  and  substitution 
of  that  function  into  Equation  (16).  For  reasons  stated  previously,  the  basic  form  of 
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the  electric-field  function  given  in  Equation  (4)  will  be  retained,  although  in  a  form 
modified  for  use  with  a  cylindrical  anode.  The  plasma  electric  field  begins  at  the 
anode  surface,  which  for  a  cylindrical  anode  is  at  r=i;,  (as  opposed  to  r=0.  the  anode 

center).  Thus,  the  presumed  electric  field  ftmction  is 


E(r)=  E.  exd 


B 

((r-r„)*b)-_ 


(17) 


where  i;,  is  the  anode  radius,  B  and  b  are  constants  that  give  specific  shape  to  the 
field,  and  the  range  of  r  is  from  i;,  to 

Following  substitution  of  E(r)  into  Equation  16.  the  resulting  lengthy  expression 
is  made  dimensionless  through  use  of  the  relations  and  definitions  given  in  Figure  M. 


,  E/E,  -;E./E,  )exp 


(fH)- 


K‘  = 


K- 


J/eD^ 


n.  = 


eE.  b 
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I'Ve' 

‘•niT 
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Figure  37  Cylindrical  Nondimensionalized  Parameters 


Note  that  q  is  the  cylindrical  counterpart  to  z,  and  that  f  has  a  more  conventional 
range  of  zero  to  infinity.  The  new  parameter  y  appears  because  the  simplifications 
that  were  previously  employed  (based  on  field  and  current  properties  at  the  electrode 
surface  and  in  the  undisturbed  plasma)  in  the  Cartesian  derivation  [Ref  2]  cannot 
eliminate  all  of  the  same  constants  when  applied  to  the  cylindrical  case,  y  has  units  of 
length,  and  can  be  thought  of  as  a  'characteristic  radius'  of  the  plasma  conditions  and 
anode  width  that  define  q  and  q  . 
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Further  simplification  of  the  differential  equation  is  achieved  with  the  derived 
expression  [Ref.  2]: 


( ^  ] 

( en  E  ' 

kX.  J 

(18) 


Finally,  order-of-magnitude  analysis  reveals  that  some  of  the  surviving  terms 
cannot  be  neglected  unless  the  anode  radius  (t;,)  is  sufficiently  large.  As  a  result,  the 

differential  equation  can  take  one  of  two  possible  forms;  the  appropriateness  of  either 
is  a  function  mainly  of  the  anode's  radial  magnitude.  Thus,  the  governing  differential 
equation  and  boundary  condition  for  'normal-sized'  anodes  (e.g.,  i;,  “10  mm)  are 

given  in  Figure  38,  while  the  corresponding  equations  for  a  'wire-thin'  anode  (e.g.. 
i;,  “0.1  mm)  are  given  in  Figure  39. 
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Figure  38  Simplified  Nondimensional  Differential  Equation 
and  Boundary  Condition,  Governing  the  1-D  'Normal  .Anode' 

Cylindrical  Problem 
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Figure  39  Simplified  Nondimensional  Differential  Equation 
and  Boundary  Condition,  (governing  the  1-D  'Wire  Anode' 

Cylindrical  Problem 
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The  equations  of  Figure  39  are  actually  appropriate  for  aU  magnitudes  of  i;,.  but 
some  of  the  terms  become  negligible  for  the  larger  anodes.  In  fact,  the  degree  of 
influence  exerted  by  the  extra  terms  of  that  differential  equation  (Figure  39)  is 
directly  tied  to  the  ratio  y/fo  ;  unless  the  cylindrical  anode's  radius  is  smaller  than,  or 

within  an  order-of-magnitude  of,  the  'characteristic  radius'  (y).  then  those  extra 
terms  are  insignificant.  For  this  reason,  the  simplified  differential  equation  of  Figure 
38  is  also  provided. 

Note  also  that  this  'normal  anode'  differential  equation  is  identical  in  form  to  the 
corresponding  planar  equation  (Figure  3).  which  indicates  that  in  spite  of  a  tortuous 
derivation,  most  cylindrical  and  planar  anodes  disturb  the  plasma  in  nearly  the  same 
manner.  Finally,  as  t,  in  the  expression  of  Figure  39  increases  toward  infinity 

(approaching  the  planar  case),  the  extra  terms  vanish  and  the  planar  differential 
equation  is  recovered. 

The  parameters  q  and  have  comparable  meaning  to  their  planar  counterparts, 
and  their  specific  values  are  computed  using  similar  techniques.  It  can  be  shown 
through  extensive  manipulations  that  the  cylindrical  equivalents  to  Equations  (5)  and 
(6)  are  as  follows: 


u  r  -'i 

b-  I  \  9  exp,2q  ) 

(19) 

Ad),  =E,VBf(q) 

(20a) 

where  f!q).q^ 

(20b) 
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Note  that  tlie  infinite  series  converges  for  all  values  of  q.  and  that  the  \/r„  term 
in  the  denominator  of  Equation  (19)  is  only  relevant  for  the  'wire-thin'  (and  smaller) 
anodes.  The  above  equations  are  combined  to  generate  the  two  important  equations 
of  Figure  40.  which  are  used  to  compute  the  parameters  q  and  q  in  accordance  with 

the  following  familiar  procedure: 

•  the  gas  composition  defines 

•  the  particular  case  or  application  specifies  E  and  n 

•  the  device  in  use  fixes  the  anode  radius  i;, 

•  q  is  then  computed  using  the  implicit  equation  at  the  top  of  Figure  40 

•  the  plasma  temperature  X,  then  yields  q  via  the  bottom  equation  of  Figure  40 

•  for  the  extremely  thin  anodes,  y  is  computed  from  its  definition  and  the  value  of 

E(f  =  0);  both  expressions  are  found  in  Figure  37. 


Figure  4(1  Equations  Thai  ^  ield  q  and  qi 


As  with  the  planar  cases,  specific  values  of  q  and  q  allo\\  computation  of  tlie 
electric  field  parameters  B  and  b  (Equations  19  and  20).  which  in  turn  yield  the  fitted 

form  of  the  E  distribution  (Equation  17).  and  ultimately  permit  the  creation  of  K\ 

/ 

(K  )  .  and  h^  profiles  in  the  plasma  regions  close  to  the  anode.  All  that  remains  to 
be  derived  are  the  cylindrical  forms  of  the  equations  that  compute  those  n-profiles, 
which  are  produced  through  manipulation  of  Equations  (2).  (13).  and  (15).  The 
resulting  elaborate  expressions  are  presented  in  Figures  41  and  42. 
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Figure  42  Electron  Population  Equation 


Gratify ingly,  these  n-profile  equations  produce  their  planar  counterparts  as  t, 
approaches  infinity.  Note  also  that  the  expressions  for  h  and  h,  are  valid  for  all 
cylindrical  anodes.  Unfortunately,  only  one  term  in  Figure  41  could  be  neglected  for 
larger,  'normal-sized'  anodes,  which  stifles  the  notion  of  a  separate  'simplified' 
equation  for  such  cases. 

B.  NUMERICALLY  SOLVED  CASES 
1.  Procedure 

The  previous  derivations  allow  detailed  numerical  analysis  of  cylindrical 
anode  sheaths  in  certain  plasma  conditions,  using  a  procedure  similar  to  that 
performed  on  the  planar  problem.  As  before,  the  governing  equation  is  rewritten  in 
terms  of  w,  as  defined  in  Equation  (21).  The  resulting  differential  equations  and 
boundary  conditions  for  both  the  'normal-sized'  anode  (~5  mm  radius  and  larger;  this 
form  also  applies  when  1)  and  the  'wire-thin'  anode  (-0.5  mm  radius  and 

smaller,  or  when  ( y/r,,)  “  1 )  cases  are  presented  in  Figures  4.1  and  44.  Note  that  the 
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form  of  the  equations  for  the  'normal-sized'  anode  is  identical  to  that  of  the 
corresponding  planar  equations  (Figure  6). 
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Fijiiire  43  Modified  1-D  Cylindrical  Equations 
for  Numerical  Solution  ('Normal-Sized'  Anode) 
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Figure  44  Modified  1-D  Cylindrical  Equations 
for  Numerical  Solution  ('Wire-Thin'  Anode) 

As  before,  application  of  a  FORTRAN  Runge-Kutta  algorithm  (see  Appendix 
A),  combined  with  the  previously-outlined  data  manipulations,  yield  the  numerical 

I 

profiles  for  K\  (K  ) ,  and  ;  for  any  desired  plasma  and  anode  width. 

2.  The  Nitrogen  Problem  with  a  Cylindrical  Anode 

The  nitrogen  plasma  conditions  designated  previously  as  'Case  F  for  the 
planar  anode  are  now'  analyzed  for  both  the  'normal-sized'  and  'wire-thin'  cylindrical 
anodes.  The  Case  I  conditions  are  repeated  on  the  next  page; 


•  nitrogen's  anode  potential  drop  (A(t), ):  15.5 1  (singly  ionized) 

•  the  electric  field  strength  in  the  undisturbed  plasma  (E. );  12.(X)0.  V/m 

•  the  particle  density  (n. ):  lO  '  m  ' 

•  the  temperature  (TJ,); 

•  the  anode  radii  under  consideration  (i;,):  10.0  mm  and  0.1mm 

These  specifications  and  the  newly-derived  equations  produce  the  following 
parameter  values  for  their  respective  anodes; 

TABLE  1  CYLINDRICAL  PARAMETERS  FOR  CASE  I  CONDITIONS 


i;,  =  10.  mm 

II 

p 

3 

3 

q 

1.75656 

1.78774 

99.1571 

102.818 

Y 

N/A 

3.9597x10' m 

b 

1.953x10"  m 

1.813x10"  m 

E 

262,543.  V/m 

293,216.  V/m 

A  comparison  reveals  that  the  parameter  values  for  the  normal-sized'  (10. 
mm)  anode  differ  negligibly  (only  0. 1  %  and  less)  from  their  corresponding  planar 
values.  This  is  not  an  unexpected  result,  in  light  of  the  cylindrical  equations' 
similarity  to  their  planar  counterparts  for  the  bigger  anodes.  More  surprising  are  the 
values  computed  for  the  'wire-thin'  (0.1  mm)  anode;  while  their  differences  from  the 
planar  values  are  not  insignificant,  only  E,  differs  by  more  than  8%.  The  extreme 

sensitivity  of  the  parameter  q,  however,  implies  that  percentage  comparisons  may  be 
misleading. 
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Direct  comparison  of  tlie  planar  and  cylindrical  profiles  best  illustrates  the 
near-negligible  differences  in  their  respective  sheath  and  ambipolar  regions.  Figures 
45  and  46  present  the  Case  I  K-tenn  profdes  produced  by  tlie  three  anode  types 
(plane,  'normal'  (10.  mm)  cylinder,  and  'wire'  (0.1  mm)  cylinder),  while  Figure  47 
depicts  the  n-profiles  for  each  of  those  anodes.  As  can  be  seen,  profiles  for  the 
'normal'  cylindrical  anode  are  visually  indistinguishable  from  the  planar  anode 
profiles.  Amazingly,  even  the  'wire'  anode  profiles  are  nearly  identical  to  those  of 
the  planar  anode.  TTiis  is  a  welcome  result;  despite  many  differences  in  the  derivation 
and  appearance  of  the  governing  equations  for  each  type  of  anode,  it  seems  that  for 
this  'small  y'  family  of  plasmas  those  cylindrical  anodes  whose  radius  is  greater  than 
that  of  a  human  hair  can  be  treated  to  behave  as  planar.  For  such  plasmas,  only  those 
cylindrical  anodes  whose  radii  are  on  the  order  of  the  sheath  thickness  itself  would 
exhibit  significant  non-planar  characteristics. 

To  further  validate  this  hypothesis,  the  three  anodes  are  also  compared  under 
tlie  lower-temperature  'Case  Ilf  circumstances,  which  represent  the  same  conditions 
as  Case  I  except  that  X.  =  300 ’K.  Consequently,  the  only  parameter  listed  in  Table  1 
that  changes  is  q  ;  its  Case  III  values  are  1983.14  for  tlie  'nonnal-sized'  cylinder,  and 

2056.36  for  the  'wire'  anode.  Figures  48.  49.  and  50  show  the  K  .  [K  ) .  and 
profiles  (respectively)  for  Case  III  conditions.  Once  again,  the  planar  and  10-mm- 
radius  anodes  disturb  the  plasma  in  identical  fashion,  while  the  0. 1-mm-radius  anode 
exhibits  nearly  negligible  variations. 

It  is  worth  noting  that  the  ratio  (v/r„)  can  theoretically  become  a  meaningful 
factor  even  for  the  'normal-sized'  anodes;  extremely  low-density  plasmas  could 
possibly  generate  values  of  y  in  the  millimeter  range.  Thus,  the  notion  of  sheath 
characteristics  being  strictly  a  function  of  anode  radius  is  valid  only  for  the 
conditions  being  considered  in  this  work. 
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Figure  45  K"  Profiles,  a  Three-Anode  Comparison 
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Figure  46  (K  )  Profiles,  a  Three-Anode  Comparison 
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Figure  47  n,,  Profiles,  a  Three-Anode  Comparison 
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Figure  48  K*  Profiles,  a  Three-Anode  Comparison 
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IV.  CONCLUSIONS 


The  manner  in  which  one-dimensional  planar  and  cylindrical  anodes  disturb 
collisional,  isothermal  plasmas  has  been  explored  in  detail.  Appropriate  sheath  and 
ambipolar  regions  have  been  shown  to  develop,  for  a  given  specific  and  reasonable 
electric  field  distribution,  across  a  wide  range  of  plasma  types  and  conditions. 
Numerical  solutions  of  the  nondimensional  governing  equations  have  allowed  specific 

and  important  observations  to  be  made  concerning: 

•  the  effect  of  temperature  variation  on  plasmas  in  the  anode  region 

•  traits  that  are  characteristic  of  various  plasma  'families',  and  traits  that  appear 

to  be  common  to  all  collisional,  low-temperature  plasmas 

•  the  effect  that  differing  anode  radii  have  on  the  sheath  for  cylindrical  anodes 

•  the  degree  of  difference  between  sheaths  of  planar  and  cylindrical  anodes 

In  addition,  an  analytical  method  has  been  presented  that  offers  a  simplified  yet 
accurate  algebraic  solution  for  the  lower-temperature  plasmas. 

The  following  briefly  summari2es  the  specifics  for  each  of  these  results. 

/ 

Appropriate  Sheath  Existence:  In  each  case  considered,  the  (K')  curv'es  showed 
that  the  net  rate-of-charge-production  decreases  and  vanishes  toward  the  outer  edge 
of  the  sheath.  Moreover,  each  of  the  n-profiles  clearly  depicted  sheath  and  ambipolar 
regions  which  were  reasonably  located  with  respect  to  the  anode  surface. 

Overall-Temperature  Variation  Effects:  Results  of  this  work  indicate  that  at 
lower  temperatures,  charge  production  in  the  outer  sheath  is  generic  to  the  electric 
field  distribution,  and  is  thus  independent  of  the  details  of  ionization  and 
recombination.  Additionally,  the  .  profiles,  and  thus  the  sheaths  themselves,  have 

been  shown  to  be  nearly  unaffected  by  substantial  changes  in  temperature.  This 
somewhat  unexpected  result,  which  held  for  every  plasma  case  considered,  reduces 
any  consequences  caused  by  the  isothennal-particle  assumption  (i.e..  the  assertion  that 
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T  =  1  =  T, ).  It  is  interesting  to  note  tliat  the  large  magnitudes  and  physical  meaning 
of  q  (the  ratio  of  electrical  energy  to  thermal  energy  in  the  sheath)  actually  predict 
such  independence  from  temperature  variations. 

Plasma  'Family'  Characteristics:  All  low-temperature,  collisional  plasmas  can  be 
partitioned  into  broad  groupings  based  on  their  specific  nondimensional  'z  or  'q' 
value  (z  denoting  a  planar  anode  problem,  and  q  signifying  cylindrical  anodes).  Each 
plasma  'family'  thus  exhibits  its  own  characteristic  electric  field  and  sheath 
properties,  even  though  its  members  may  be  widely  diverse  in  composition  and 
condition. 

Plasmas  with  large  z  (or  q)  values  show  significant  charge  production  throughout 
the  entire  sheath,  in  contrast  to  the  localized  and  decreased  production  rates  depicted 
for  the  'small  z/q'  plasmas.  This  may  indicate  the  fact  that  electric-field  effects  are 
more  modest  and  thus  localized  in  those  plasmas  with  smaller  values  of  z  or  q.  In 
addition,  the  magnitude  of  the  charged  particles  near  the  anode  changes  dramatically 
as  'z/q'  varies,  especially  in  the  sheath.  The  sheaths  of  'large  z/q'  plasmas  exhibit 
charged-particle  densities  up  to  one  order-of-magnitude  smaller  than  those  of  the 
'small  z/q'  sheaths.  More  importantly,  the  sheath  itself  extends  further  from  the 
anode  surface  as  the  value  of  'z/q'  decreases.  In  particular,  the  sheath  of  the  'small 
z/q'  plasmas  stretches  over  ten  times  the  distance  occupied  by  the  'large  z/q'  sheath. 

Cylindrical  Anode  Radius  Effects:  The  sheatii  and  the  charge  production  rates 
have  been  shown  to  vary  as  the  anode  radius  varies;  however,  the  differences  are 
almost  negligible  for  all  practical  anodes  and  conditions.  It  appears  that  only 
extremely  low-density  plasmas,  or  cylindrical  anodes  whose  radii  are  on  the  order  of 
the  sheath  thickness  itself  would  exhibit  significantly  different  characteristics. 
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Planar  and  Cylindrical  Anode  Differences:  In  an  extension  of  the  previous 
paragraph,  the  results  of  this  work  indicate  that  for  the  range  of  conditions  explored, 
all  practical  cylindrical  anodes  can  be  treated  to  behave  as  planar  anodes. 

Simple  Analytical  Approximations  fLow'  Temperatures  Only):  Using  the  'Outer 
Solution'  technique,  an  algebraic  solution  can  be  used  to  generate  highly  accurate 
approximations  to  low-temperature  anode  problems.  In  addition,  as  long  as 
derivatives  of  the  presumed  electric  field  distribution  exist,  the  sheath  profiles  for  all 
cases  can  be  produced  analytically  to  any  desired  accuracy. 

One  suggestion  for  further  work  involves  researching  the  ef.^cts  that  other 
electric  field  distributions  w'ould  have  on  the  anode  problem.  Such  distributions  may 
be  derived  from  empirical  data  or  from  other  'guesstimates',  and  may  have  three  or 
more  adjustable  parameters.  The  most  prominent  reason  for  this  suggestion  is  the 
need  to  corroborate  the  generic  charge  production  rate  observed  at  lower 
temperatures. 

Another  area  requiring  further  study  is  the  exploration  of  two-temperature 
plasmas  and  their  anode  sheaths.  Results  for  such  conditions  could  confirm  the 
assertion  of  sheath  insensitivity  to  temperature  changes. 
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APPENDIX  A 


Applicable  FORTRAN  Programs 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3.1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


DOUBLE  PRECISION  YA ( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 
DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W, KPLUS , WPRM , KPLSPRM 
REAL  PI ,XL,XCHK,EINF,EO,Z,ETA 
INTEGER  NS, PI 
PRINT  * 

PRINT  FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 

PRINT  FOR  THE  CARTESIAN  COMPUTATION  ' 

PRINT  OF  W,  KPLUS,  S,  KPLSPRM' 

PRINT  *, 'INPUT  VALUES  FOR  EINF,  EO,  Z,  AND  ETA' 

READ  *, EINF, E0,Z, ETA 

1  IM=1  !  Number  of  equations 

Y(l)  =  1.00  I  Initial  condition  for  y~l  at  x=XP. 

C  Y(2)  =  0  !  Initial  condition  for  y~2  at  x=XP  (if  nec) 

PRINT  * 

PRINT  *, 'INTERVAL  OF  X  FOR  PRINTING  ?' 

READ  *,PI 

PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-x)' 

READ  *,H 

NS  =  NINT(PI/H) 

PRINT  *, 'MAXIMUM  X  TO  STOP  CALCULATION  ?' 

READ  *,XL 

PRINT  * , '  H=  ' ,H 

Pl  =  0 
XP=0 
HH=H/2 
KPLUS=1 . 0 


print  *,'  HH=  ' , HH 
print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0 

PRINT* , ' LINE  y/a 
WRITE  (*,98)  LI,XP 

28  LI=LI+1 

DO  N=1,NS 
XB  =  XP 
XP=XP+H 
XM=XB+HH 

J=1 

DO  1=1, IM 

YA(I)=Y(I) 
END  DO 
XA=XB 


!  Line  no.  initialization 
KPLUS' 

KPLUS 


!  Old  time 
1  New  time 
I  Midpoint  time 

!  This  part  computes  k~l. 
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CALL  FUNCT( EK, J, YA,H,XA,EINF,EO,Z,ETA) 

J=2  !  This  part  computes  k~2. 

DO  1=1, I M 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT(EK, J, YA,H,XA, EINF,E0,Z,ETA) 

J=3  !  This  part  computes  k~3. 

DO  1=1, IM 

YA(I)=Y(I)+EK(2,I)/2 
END  DO 
XA=XM 

CALL  FUNCT(EK, J, YA,H,XA,EINF,EO,Z,ETA) 

J=4  !  This  part  computes  k''4. 

DO  1=1, IM 

YA( I )=Y( I )+EK(  3,1) 

END  DO 
XA=XP 

CALL  FUNCT(EK, J, YA,H,XA,EINF,EO,Z,ETA) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y( I ) =Y( I ) + ( EK ( 1 , I ) +EK (2,1 ) *2+EK (3,1 ) *2+EK ( 4 , I ) ) /6 
END  DO 


C 

c 


c 

c 


C 

c 

c 


W=Y( 1 ) 

You  now  have  W(  y/a  )  --  to  get  K+(  y/a  )  you  must 
multiply  by  {  E(  y/a  )  /  E(0)  } 

E=( EINF/EO ) *DEXP( (Z**2.  )/( (XP+1. )**2. ) ) 

KPLUS=W*E 

To  get  the  derivative  of  W  [  W' (  y/a  )  ),  just  plug  the 

computed  values  of  W  back  in  the  original  Ist-order  ODE 


J=5 

YA ( 1 ) =W 
XA=XP 

CALL  FUNCT(EK, J,YA,H,XA,EINF,EO,Z,ETA) 
WPRM=EK( 5,1 )/H 


NOW,  to  get  d/dy  ( K+ 1  from  d/dy  I  K+/(E/E0)  ]  ,  must 

perform  the  following  operation: 

KPLSPRM=(WPRM*(E**2. ) +KPLUS*E* ( -2 . * ( Z**2 . )/( (XP+1. )**3. ) ) )/E 


To  keep  from  generating  unplottable  50,000  point  data  files, 
the  following  will  edit  out  data  points  depending  on 
'where'  they  occur 


XCHK=(XP/H) 

IF  (XCHK  .GT. 
IF  (XCHK  .GT. 
IF  (XCHK  .GT. 
IF  (XCHK  .GT. 
GOTO  4  4 


100000)  GOTO  72 
10000)  GOTO  73 
1000)  GOTO  74 
100)  GOTO  75 


72 


Pl=Pl+l 


72 


73 

74 

75 


IF  (PI  .NE. 
Pl  =  0 
GOTO  44 
Pl=Pl+l 
IF  (PI  .NE. 
P1  =  0 
GOTO  4  4 
Pl=Pl+l 
IF  (PI  .NE. 
P1  =  0 
GOTO  44 
P1=P1+1 
IF  (PI  .NE. 
P1  =  0 


10000)  GOTO  88 


1000)  GOTO  88 


100)  GOTO  88 


10)  GOTO  88 


44  WRITE  (13,*)  XP,  KPLSPRM 

WRITE  (12,*)  XP,  KPLUS 
C  WRITE  (10,*)  XP,  W 

88  CONTINUE 


END  DO 

WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  I2,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT.  XL)  GOTO  28 

200  PRINT* 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP.' 
READ  * , K 

IF(K.EQ.l)  GOTO  1 
PRINT* 

END 


SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA )  !  DEFINES  SET  OF  EQS 

DOUBLE  PRECISION  EK ( 0 : 4 , 0 : 1 0 ) , YA( 0 : 1 0 ) , H , XA , PARTl ( 0 : 4 , 0 : 1 0 ) 

DOUBLE  PRECISION  PART2 ( 0 : 4 , 0 : 1 0  ) 

PARTI ( J,1 )=DEXP( (2.*Z**2.)/( (XA+1. )**2. ) )/( (XA+1. )**3. ) 

PART2( J , 1 )=( DEXP( ( Z**2 . )/( (XA+1 . ) **2 . ) ) ) * ( ETA*EINF/E0 ) 

EK( J,1)=(2.*ETA-PART2( J,1 )*YA( 1 ) - ( ETA* ( EINF/EO ) **2 . )*PARTl( J,l) ) *H 
C  EK(J,2)=  the  second  ode,  if  nec. 

RETURN 

END 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3.1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


DOUBLE  PRECISION  YA( 0 : 1 0 ) , YN { 0 : 10 ) , EK ( 0 : 4 , 0 : 10 ) , Y ( 0 ; 10 ) , XA 

DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 

DOUBLE  PRECISION  NI , NE , NI 1 , NI 2 , NI 3 , NI2A , NI 3A, NEl 

REAL  PI ,XL,XCHK,EINF,EO,Z,ETA 

REAL  EPSO,CHRG,K,TEMP,NINF,A 

INTEGER  NS, Pi 

EPS0=8.853742E-12 

CHRG=1.602E-19 

K=1.38054E-23 

PRINT  * 

PRINT  * , '  FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 

PRINT  FOR  THE  CARTESIAN  COMPUTATION  ' 

PRINT  OF  NI  AND  NE ' 

PRINT  *, 'INPUT  VALUES  FOR  EINF,  EO,  Z,  AND  ETA' 

READ  *, EINF, E0,Z, ETA 

PRINT  *, 'INPUT  VALUES  FOR  TEMP.  AND  NINF' 

READ  *, TEMP, NINF 

PRINT  *, 'INPUT  VALUE  FOR  a' 

READ  *,A 

1  IM=1  !  Number  of  equations 

Y(l)  =  1.00  !  Initial  condition  for  y~l  at  x=XP. 

C  Y(2)  =  0  !  Initial  condition  for  y~2  at  x=XP  (if  nec) 

PRINT  * 

PRINT  *, 'INTERVAL  OF  X  FOR  PRINTING  ?' 

READ  *,PI 

PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-x)' 

READ  *,H 

NS  =  NINT(PI/H) 

PRINT  *, 'MAXIMUM  X  TO  STOP  CALCULATION  ?' 

READ  *,XL 

PRINT  * , '  H=  ' ,H 

P1  =  0 
XP=0 
HH=H/2 
KPLUS=1 . 0 


print  * , '  HH=  ' ,HH 
print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0  !  Line  no.  initialization 

PRINT*, 'LINE  y/a  KPLUS' 

WRITE  (*,98)  LI,XP,  KPLUS 
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28 


LI=LI+1 
DO  N=1 ,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 


!  Old  time 
I  New  time 
!  Midpoint  time 


J=1  !  This  part  computes  k~l. 

DO  1=1, IM 

YA(I)=Y(I) 

END  DO 
XA=XB 

CALL  FUNCT( EK, J, YA,H,XA,EINF,EO,Z,ETA) 

J=2  !  This  part  computes  k~2. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA ) 

J=3  !  This  part  computes  k~3. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 2 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT(EK, J,YA,H,XA, EINF, E0,Z, ETA) 

J=4  !  This  part  computes  k~4. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 3 , I  ) 

END  DO 
XA=XP 

CALL  FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA ) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2 , I )*2+EK( 3, I ) *2+EK( 4 , I ) )/6 
END  DO 


W=Y( 1 ) 

c  You  now  have  W(  y/a  )  --  to  get  K+(  y/a  )  you  must 

c  multiply  by  {  E(  y/a  )  /  E(0)  } 

E=( EINF/EO ) *DEXP(  (Z**2.  )/(  (XP+1.  )**2.  )  ) 

KPLUS=W*E 

c  To  get  the  derivative  of  W  [  W'(  y/a  )  ],  just  plug  the 

c  computed  values  of  W  back  in  the  original  Ist-order  ODE 

J=5 

YA( 1 )=W 
XA=XP 

CALL  FUNCT(EK, J,YA,H,XA, EINF, E0,Z, ETA) 

WPRM=EK( 5 , 1 )/H 

NOW,  to  get  d/dy  [ K+ ]  from  d/dy  [  K+/(E/E0)  ]  ,  must 

perform  the  following  operation: 

KPLSPRM=( WPRM* ( E**2 . )+KPLUS*E* (-2.*(Z**2.)/(( XP+1 , ) **3 . ) ) )/E 


75 


n  n  o  o  n  n  o 


AND  NOW,  to  generate  the  non-dimensionalized  n-curves 
ne(  y/a  )  and  ni(  y/a  ) 
just  compute.'  these  equations: 

NI1= ( KPLUS/E ) * ( EINF/EO ) 

NI2A=( EPSO*K*TEMP)/( (CHRG**2. ) *NINF* ( A* * 2 . ) ) 

NI2=NI2A* ((6.*(Z**2.)/( (XP+1 .)>*4.))+(4.*(Z**4.)/( (XP+1 . ) **6 . ) ) ) 

NI3A=( ( EPS0*E0 )/( CHRG*F :NF*A) ) *E 
NI3=NI3A* ( -2 . * ( Z**2 . )/\ ;XP+1. )**3. ) ) 

NI=.5*(NI1+NI2+NI3) 

NE1=(  (E0*EPS0*(Z**2.  )  )/(CHRG*NINF*A) )*E*(2./(  (XP+1.  )**3.  )  ) 
NE=NI+NE1 


To  prevent  every  data  point  from  being  written  to  the  file 
(resulting  in  unplottable  50,000  pt.  data  files),  the  following 
edits  out  a  percentage  of  the  data  depending  on  'where'  it 
was  generated. 

XCHK=(XP/H) 

IF  (XCHK  .GT.  100000)  GOTO  72 
IF  (XCHK  .GT.  10000)  GOTO  73 
IF  (XCHK  .GT.  1000)  GOTO  74 
IF  (XCHK  .GT.  100)  GOTO  75 
GOTO  44 

72  Pl=Pl+l 

IF  (Pi  .NE.  10000)  GOTO  88 
Pl  =  0 
GOTO  44 

73  P1=P1+1 

IF  (PI  .NE.  1000)  GOTO  88 
P1  =  0 
GOTO  44 

74  P1=P1+1 

IF  (PI  .NE.  100)  GOTO  88 
Pl  =  0 
GOTO  44 

75  P1=P1+1 

IF  (Pi  .NE.  10)  GOTO  88 
P1  =  0 


44  WRITE  (10,*)  XP,  NI 

WRITE  (11,*)  XP,  NE 

88  CONTINUE 


END  DO 

WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  12,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT.  XL'  GOTO  28 

200  PRINT-' 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP.' 
READ  *,K1 

IF(Kl.EQ.l)  GOTO  1 
PRINT* 
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SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , EO , Z , ETA )  !  DEFINES  SET  OF  EQS 

DOUBLE  PRECISION  EK { 0 : 4 , 0 : 1 0 ) , YA( 0 : 1 0 ) , H , XA , PARTI ( 0 : 4 , 0 : 1 0 ) 

DOUBLE  PRECISION  PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTI ( J, 1 )=DEXP( (2.*Z**2.)/(( XA+1 . ) **2 . ) )/( (XA+l . ) **3 . ) 

PART2 ( J , 1 ) = ( DEXP( ( Z**2 . )/( ( XA+1 . ) **2 . ) ) ) * ( ETA*EINF/EO ) 

EK( J, 1 )=( 2 . *ETA-PART2 ( J , 1 ) *YA( 1 )-( ETA* ( EINF/EO ) **2 . ) * PARTI (J , 1 ) ) *H 
EK(J,2)=  the  second  ode,  if  nec. 

RETURN 

END 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3.1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


1 

C 


DOUBLE  PRECISION  YA( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , 
DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 
REAL  PI,XL,XCHK,EINF,ERO,Q,ETA 
INTEGER  NS, PI 
PRINT  * 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 

PRINT  *  ,  ' 


FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 
FOR  THE  CYLINDRICAL  COMPUTATION  ' 

OF  W,  KPLUS,  &  KPLSPRM' 

(  NORMAL-SIZED  ANODE:  rO  =  1  cm  ) ' 


PRINT  *, 'INPUT  VALUES  FOR  EINF,  ErO,  Q,  AND  ETA' 

READ  *, EINF, ERO, Q, ETA 

IM=1  !  Number  of  equations 

Y(l)  =  1.00  !  Initial  condition  for  wl  at  r~=XP. 

Y(2)  =  0  1  Initial  condition  for  w2  at  r~=XP  (if  nec) 


PRINT  * 

PRINT  *, 'INTERVAL  OF  r~  FOR  PRINTING  ?' 
READ  *,PI 

PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-r~)' 
READ  *,H 


NS  =  NINT(PI/H) 

PRINT  *, 'MAXIMUM  r~  TO  STOP  CALCULATION  ?' 
READ  *,XL 

PRINT  *  ,  '  dr*"  =  '  ,H 


P1  =  0 
XP=0 
HH=H/2 , 
KPLUS=1 . 0 


print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0  !  Line  no.  initialization 

PRINT*, 'LINE  r~  KPLUS' 

WRITE  (*,98)  LI,XP,  KPLUS 


28 


C 


LI=LI+1 
DO  N=1 ,NS 
XB=XP 
XP=XP+H 
XM=XB+HH 


Old  time 
New  time 
Midpoint  time 


-Runge-Kutta  Scheme- 


J  =  1 

DO  1=1, I M 


This  is  part  1. 
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YA( I )=Y( I ) 

END  DO 
XA=XB 

CALL  FUNCT(EK, J, YA,H,XA,EINF,ERO,Q,ETA) 

J=2  !  This  is  part  2. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT( EK , J , YA , H , XA , EINF , ERO , Q , ETA ) 

J=3  !  This  is  part  3. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 2, I )/2 
END  DO 
XA=XM 

CALL  FUNCT( EK,J,YA,H, XA, EINF, ERO, Q, ETA) 

J=4  !  This  is  part  4. 

DO  1=1, IM 

YA( I )=y( I )+EK( 3,1) 

END  DO 
XA=XP 

CALL  FUNCT ( EK , J , YA , H , XA, EINF , ERO , Q , ETA ) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2 , I ) *2+EK ( 3 , I ) * 2+EK ( 4 , I ) )/6 
END  DO 


W=Y( 1 ) 

You  now  have  W(  r~  )  —  to  get  K+(  r~  )  you  must 
multiply  by  {  E(  r~  )  /  E(0)  ) 

E=( EINF/ERO ) *DEXP( (Q**2 .  )/(  (XP  +  1 .  )**2  .  ) ) 

KPLUS=W*E 

To  get  the  derivative  of  W  (  W'(  r~  )  ],  just  plug  the 

computed  values  of  W  back  in  the  original  Ist-order  ODE 

J=5 

YA( 1 )=W 
XA=XP 

CALL  FUNCT(EK, J,YA,H,XA, EINF, ERO ,Q, ETA) 

WPRM=EK(5,1)/H 

NOW,  to  get  d/dr~  [ K+ ]  from  d/dr~  (  K+/(E/ErO)  ]  ,  must 

perform  the  following  operation: 

KPLSPRM=( WPRM* ( E**2 .  )+KPLUS*E* (-2.*(Q**2.)/(( XP+1 .  )  *  *  3  .  )  )  ) /E 

To  keep  from  generating  unplottable  50,000  point  data  files, 
the  following  edits  out  some  of  the  data  points  depending 
on  'where'  they  occur 


XCHK=(XP/H) 

IF  (XCHK  .GT.  100000)  GOTO  72 
IF  (XCHK  .GT.  10000)  GOTO  73 
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IF  (XCHK  .GT.  1000)  GOTO  74 
IF  (XCHK  ,GT.  100)  GOTO  75 
GOTO  44 


72 

Pl=Pl+l 

IF  (Pi  .NE. 
Pl  =  0 

GOTO  44 

10000 

)  GOTO  88 

73 

Pl=Pl+l 

IF  (PI  .NE. 
P1  =  0 

GOTO  44 

1000) 

GOTO  88 

74 

Pl=Pl+l 

IF  (PI  .NE. 
P1  =  0 

GOTO  44 

100) 

GOTO  88 

75 

P1=P1+1 

IF  (PI  .NE. 

10)  GOTO  88 

P1  =  0 

44 

WRITE  (13, 

*)  XP, 

KPLSPRM 

WRITE  (12, 

*)  XP, 

KPLUS 

WRITE  (14, 

*)  XP, 

W 

88  CONTINUE 


END  DO 

WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  I2,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT,  XL)  GOTO  28 

200  PRINT* 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP.' 

READ  *,K1 

IF(Kl.EQ.l)  GOTO  1 
PRINT* 

END 

SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA )  !  DEFINES  SET  OF  EQS 

DOUBLE  PRECISION  EK ( 0 : 4 , 0 : 1 0 ) , YA( 0 : 10 ) , H , XA , PARTI ( 0 : 4 , 0 : 1 0 ) 

DOUBLE  PRECISION  PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTI ( J , 1 )=ETA* ( DEXP( ( 2 . *Q**2 . )/( (XA+1 . ) **2 . ) )/( (XA+1 . ) **3 , ) ) 
PART2( J,1)=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) ) * ( ETA*EINF/ER0 ) 

EK( J,1)=(2.*ETA-PART2( J,l)*YA(l)-( ( EINF/ERO ) * *2 . )*PARTl( J,l) ) *H 
C  EKd,2)=  the  second  ode,  if  nec. 

RETURN 

END 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3,1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


DOUBLE  PRECISION  YAl 0 : 10 ) , YN ( 0 : 1 0 ) , EK ( 0 ; 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) , XA 

DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 

DOUBLE  PRECISION  NI , NE , NI 1 , NI 2 , NI 3 , NEl 

REAL  CC1,CC2,CC3,S2,S3,S2A,S2B,S2C 

REAL  PI ,XL,XCHK,EINF,ER0,Q,ETA 

REAL  EPSO ,CHRG,K,TEMP,NINF,B,R0 

INTEGER  NS, PI 

EPS0=8 . 853742E-12 
CHRG=1 .602E-19 
K=1.38054E-23 
R0=.01 

PRINT  * 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 


FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 
FOR  THE  CYLINDRICAL  COMPUTATION  ' 

OF  NE  AND  NI' 

(  NORMAL-SIZED  ANODE:  rO  =  1  cm  ) ' 


PRINT  *, 'INPUT  VALUES  FOR  EINF,  ErO,  Q,  AND  ETA' 

READ  *, EINF, ERO, Q, ETA 

PRINT  *, 'INPUT  VALUES  FOR  TEMP.  AND  NINF' 

READ  *, TEMP, NINF 

PRINT  *, 'INPUT  VALUE  FOR  b' 

READ  *,B 

IM=1  !  Number  of  equations 

Y(l)  =  1.00  !  Initial  condition  for  wl  at  r''“XP. 

Y(2)  =  0  I  Initial  condition  for  w2  at  r~=XP  (if  nec) 

PRINT  * 

PRINT  *, 'INTERVAL  OF  r~  FOR  PRINTING  ?' 

READ  *,PI 

PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-r~)' 

READ  *,H 


NS  =  NINT(PI/H) 

PRINT  *, 'MAXIMUM  r~  TO  STOP  CALCULATION  ?' 
READ  *,XL 

PRINT  * , '  dr~  =  ' ,H 

Pl  =  0 
XP  =  0 
HH=H/2 . 

KPLUS=1 . 

print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0  !  Line  no.  initialization 

PRINT*, 'LINE  r~  KPLUS' 

WRITE  (*,98)  LI,XP,  KPLUS 
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28 


C 


LI=LI+1 
DO  N=1,NS 
XB=XP 
XP=XP+H 
XM*XB+HH 


!  Old  time 
!  New  time 
1  Midpoint  time 


Runge-Kutta  Scheme 


J=1  !  This  is  part  1. 

DO  1=1, IM 

YA( I )=Y( I ) 

END  DO 
XA=XB 

CALL  FUNCT (EK,J,YA,H,XA,EINF, ERO , Q , ETA ) 


J=2  !  This  is  part  2. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT( EK, J, YA,H,XA,EINF,ERO,Q,ETA) 

J=3  !  This  is  part  3. 

DO  1=1, IM 

YA(I)=Y(I)+EK(2,I)/2 
END  DO 
XA=XM 

CALL  FUNCT( EK, J, YA,H,XA,EINF, ER0,Q,ETA) 


J=4  !  This  is  part  4. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 3,1) 

END  DO 
XA=XP 

CALL  FUNCT(EK, J,YA,H,XA,EINF,ERO,Q,ETA) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2, I )*2+EK( 3, I ) *2+EK( 4 , I ) )/6 
END  DO 


C 


W=Y( 1 ) 

c  You  now  have  W(  r~  )  —  to  get  K+(  r~  )  you  must 

c  multiply  by  {  E(  r~  )  /  E(0)  } 

E=( EINF/ERO ) *DEXP(  (Q**2.  )/(  (XP+1.  )**2.  )  ) 

KPLUS=W*E 

c  To  get  the  derivative  of  W  (  W'(  r~  )  Jr  just  plug  the 

c  computed  values  of  W  back  in  the  original  Ist-order  ODE 


J=5 

YA( 1 )=W 
XA=XP 

CALL  FUNCT ( EK , J , YA , H , XA , El NF , ERO , Q , ETA ) 
WPRM=EK( 5, 1 )/H  o. 


■  nnn  no 


NOW,  to  get  d/dr~  [ K+ ]  from  d/dr~  (  K+/(E/ErO)  ]  ,  must 
perform  the  following  operation: 

KPLSPRM=(WPRM* ( E**2 .  ) +KPLUS*E* ( -2 . * ( Q* *2 .  )/(  (XP+1.  ) *  *  3  .  )  )  ) /E 


AND  NOW,  to  generate  the  nondimensional i zed  n-curves 
ne(  r~  )  and  ni(  r~  ) 
compute  the  following  equations 

NI 1= ( KPLUS/E ) * ( EINF/ERO ) 

CC2= ( K*TEMP*EPSO )/( ( CHRG**2 . ) *NINF* ( B**2 . ) ) 
S2A=(6.*(Q**2. ) )/( (XP+1. )**4. ) 

S2B=(4.*(Q**4. ) )/( (XP+1.  )**6.  ) 

S2C=(1/(XP+(R0/B) ) )*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) ) 
S2=S2A+S2B+S2C 
NI2=CC2*S2 

CC3=( EPS0*ER0 )/( CHRG*N1NF*B) 

S3=(  (-2.*(Q**2.  )  )/(  (XP  +  1.  )**3.  ) )  +  (l./(XP+(R0/B)  )  ) 
NI3=CC3*E*S3 

NI=. 5* (NI1+NI2+NI3 ) 

CC1= ( EPS0*ER0 )/( CHRG*NINF*B ) 

NE1=CC1*E*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) + ( 1 ./( XP+ ( RO/B ) ) ) ) 
NE=NI-NE1 


To  keep  from  generating  unplottable  50,000  point  data  files, 
the  following  edits  out  some  of  the  data  points  depending 
on  'where'  they  occur 

XCHK=(XP/H) 

IF  (XCHK  .GT.  100000)  GOTO  72 
IF  (XCHK  .GT.  10000)  GOTO  73 
IF  (XCHK  .GT.  1000)  GOTO  74 
IF  (XCHK  .GT.  100)  GOTO  75 
GOTO  44 

72  P1=P1+1 

IF  (PI  .NE.  10000)  GOTO  88 
P1  =  0 
GOTO  44 

73  P1=P1+1 

IF  (Pi  .NE.  1000)  GOTO  88 
P1  =  0 
GOTO  44 

74  P1=P1+1 

IF  (Pi  .NE.  100)  GOTO  88 
P1  =  0 
GOTO  44 

75  P1=P1+1 

.  IF  (Pi  .NE.  10)  GOTO  88 

P1  =  0 

44  WRITE  (10,*)  XP,  NI 

WRITE  (11,*)  XP,  NE 


88 


CONTINUE 
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END  DO 


WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  12,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT.  XL)  GOTO  28 

200  PRINT* 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP.' 

READ  *,K1 

IF(Kl.EQ.l)  GOTO  1 
PRINT* 

END 

C***************************************** 

SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA )  !  DEFINES  SET  OF  EQS 

DOUBLE  PRECISION  EK ( 0 : 4 , 0 : 1 0 ) , YA( 0 : 1 0 ) , H , XA , PARTI ( 0 : 4 , 0 : 1 0 ) 

DOUBLE  PRECISION  PART2 ( 0 : 4 , 0 : 1 0 ) 

PARTI (J,l )=ETA*(DEXP( (2,*Q**2. )/( (XA+1. )**2. ) )/{ (XA+1, )**3. ) ) 
PART2( J,1)=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) ) * ( ETA* EINF/ERO ) 

EK( J,1)=(2.*ETA-PART2( J,l)*YA(l)-( ( EINF/ERO )** 2 . )*PARTl(J,l) ) *H 
C  EK(J,2)=  the  second  ode,  if  nec. 

RETURN 

END 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3.1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


DOUBLE  PRECISION  YA ( 0 : 1 0 ) , YN ( 0 : 1 0 ) , EK ( 0 : 4 , 0 : 1 0 ) , Y ( 0 : 1 0 ) 
DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 
REAL  PI ,XL,XCHK,EINF,ERO,Q,ETA,RO,G,B,EPSO,CHRG 
INTEGER  NS, PI 
PRINT  * 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 


FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 
FOR  THE  CYLINDRICAL  COMPUTATION  ' 
OF  W,  KPLUS,  &  KPLSPRM' 

(  WIRE-THIN  ANODE:  rO  =  0.1  mm  )' 


XA 


PRINT  *, 'INPUT  VALUES  FOR  EINF,  ErO,  Q,  AND  ETA' 
READ  *, EINF, ERO, Q, ETA 

PRINT  *, 'INPUT  VALUES  FOR  GAMMA  AND  B' 

READ  *,G,B 


R0=.0001 

EPS0=8.853742E-12 
CHRG=1 .602E-19 

IM=1  !  Number  of  equations 

Y(l)  =  1.00  !  Initial  condition  for  wl  at  r~=XP. 

Y{2)  =  0  !  Initial  condition  for  w2  at  r~=XP  (if  nec) 

PRINT  * 

PRINT  *, 'INTERVAL  OF  r~  FOR  PRINTING  ?' 

READ  *,PI 

PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-r")' 

READ  *,H 

NS  =  NINT(PI/H) 

PRINT  *, 'MAXIMUM  r~  TO  STOP  CALCULATION  ?' 

READ  *,XL 

PRINT  * ,  '  dr~  =  '  ,H 

Pl  =  0 
XP  =  0 
HH=H/2 . 

KPLUS  =  1  . 

print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0  !  Line  no.  initialization 

PRINT*, 'LINE  r~  KPLUS' 

WRITE  (*,98)  LI,XP,  KPLUS 


LI«LI+1 
DO  N=1,NS 
XB  =  XP 
XP=XP+H 


1  Old  time 
!  New  time 
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o  n 


XM=XB+HH 


I  Midpoint  time 


C  - Runge-Kutta  Scheme - 

J=1  I  This  is  part  1. 


DO  1  =  1,  IM 

YA( I )=Y( I  ) 

END  DO 
XA=XB 

CALL  FUNCT (EK,J,YA,H,XA,EINF, ERO , Q , ETA , G , B ) 

J=2  !  This  is  part  2. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT( EK, J, YA,H,XA,EINF,ERO,Q,ETA,G,B) 

J=3  !  This  is  part  3. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 2 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT(EK, J, YA,H,XA,EINF,ERO,Q,ETA,G,B) 

J=4  !  This  is  part  4. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 3,1) 

END  DO 
XA=XP 

CALL  FUNCT( EK, J, YA,H,XA,EINF,ERO,Q,ETA,G,B) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y( I )=Y( I )+( EK( 1 , I )+EK( 2, I )*2+EK( 3, I ) *2+EK( 4 , I ) )/6 
END  DO 


C 


W=Y( 1  ) 

c  You  now  have  W(  r~  )  --  to  get  K+(  r"  )  you  must 

c  multiply  by  {  E(  r~  )  /  E(0)  } 

E=( EINF/ERO ) *DEXP(  (Q**2.  )/(  (XP+1.  )**2.  )  ) 

KPLUS=W*E 

c  To  get  the  derivative  of  W  [  W' (  )  1 »  just  plug  the 

c  computed  values  of  W  back  in  the  original  Ist-order  ODE 


J=5 

YA ( 1 ) =W 
XA=XP 

CALL  FUNCT(EK, J, YA,H,XA,EINF,ERO,Q,ETA,G,B) 

WPRM=EK( 5 , 1 )/H 

NOW,  to  get  d/dr~  ( K+ ]  from  d/dr"  (  K+/(E/ErO)  ]  ,  must 

perform  the  following  operation: 

KPLSPRM=(WPRM* ( E**2 . ) +KPLUS*E* ( -2 . *( 0**2 . )/( ( XP+1 . ) **3 . ) ) )/E 
C  To  keep  from  generating  unplottable  50,000  point  data  files. 
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o  n 


72 


73 


74 


75 


the  following  edits  out  some  of  the  data  points  depending 
on  'where'  they  occur 


XCHK=(XP/H) 

IF  (XCHK  .GT.  100000)  GOTO  72 
IF  (XCHK  .GT.  10000)  GOTO  73 
IF  (XCHK  .GT.  1000)  GOTO  74 
IF  (XCHK  .GT.  100)  GOTO  75 
GOTO  44 


P1=P1+1 
IF  (PI  .NE. 
P1  =  0 
GOTO  44 
P1=P1+1 
IF  (PI  .NE. 
P1  =  0 
GOTO  44 
Pl=Pl+l 
IF  (PI  .NE. 
Pl  =  0 
GOTO  44 
P1=P1+1 
IF  (PI  .NE. 
P1  =  0 


10000)  GOTO  88 


1000)  GOTO  88 


100)  GOTO  88 


10)  GOTO  88 


44 

WRITE 

(13,*) 

XP, 

KPLSPRM 

WRITE 

(12,*) 

XP, 

KPLUS 

WRITE 

(14,*) 

XP, 

W 

88  CONTINUE 


END  DO 

WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  12,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT.  XL)  GOTO  28 

200  PRINT* 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP.' 
READ  *,K1 

IF(Kl.EQ.l)  GOTO  1 
PRINT* 

END 


SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B )  !  DEFINES  ODES 

DOUBLE  PRECISION  EK ( 0 : 4 , 0 : 1 0 ) , YA ( 0 ; 1 0 ) , PART2 ( 0 : 4 , 0 : 1 0 ) 

DOUBLE  PRECISION  PART3 ( 0 : 4 , 0 ; 1 0 ) , PART4 ( 0 : 4 , 0 : 1 0 ) , H , XA , Cl 

REAL  RO 

R0=.0001 

C2=ETA*G*( (EINF/ER0)**2.  ) 

PART4 ( J , 1 )=C2* ( DEXP(  (2.*(Q**2.))/( (XA+1 .  ) **2 .  )  )  )/(XA*B  +  R0 ) 

C1  =  ETA*( (G/R0)  +  1. )*(  (EINF/ER0)**2.  ) 

PART3( J,1 )=C1*(DEXP( (2.*Q**2. )/( (XA+1. )**2. ) )/( (XA+1. )**3. ) ) 
PART2( J,1)=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) ) * ( ETA*EINF/ER0 ) 

EK ( J , 1 ) = ( 2 . *ETA-PART2 ( J , 1 ) *YA( 1 ) -PART3 ( J , 1 ) +PART4 ( J , 1 ) ) *H 
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C  EK(J,2)=  the  second  ode, 

RETURN 
END 


f  nec . 
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THIS  PROGRAM  IS  BUILT  AROUND  A  GRAFkit  3.1  RUNGE-KUTTA 
ALGORITHM  PRODUCED  BY  SCO,  INC. 


1 


DOUBLE  PRECISION  YA(  0  : 1 0  )  ,  YN (  0  : 1 0  )  ,  EK (  0  :  4 , 0  ;  1 0  )  • '^  (  0  : 1 0  )  ,  XA 

DOUBLE  PRECISION  XP , XB , XM , H , HH , E , W , KPLUS , WPRM , KPLSPRM 

DOUBLE  PRECISION  NI , NE , NI 1 , NI 2 , NI 3 , NEl 

REAL  CC1,CC2,CC3,S2,S3,S2A,S2B,S2C,S2D 

REAL  PI,XL,XCHK,EINF,ER0,Q,ETA 

REAL  EPSO , CHRG , K , TEMP , NINF , B , G , RO 

INTEGER  NS, PI 

EPS0=8 . 853742E-12 
CHRG=1 .602E-19 
K=1.38054e-23 
R0=.0001 

PRINT  * 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 

PRINT  * , ' 


FOURTH  ORDER  RUNGE-KUTTA  SCHEME  ' 
FOR  THE  CYLINDRICAL  COMPUTATION  ' 
OF  NE  AND  NI ' 

(  WIRE-THIN  ANODE:  rO  =  0.1  mm  )' 


PRINT  *, 'INPUT  VALUES  FOR  EINF,  ErO,  Q,  AND  ETA' 

READ  *, EINF, ERO ,Q, ETA 

PRINT  *, 'INPUT  VALUES  FOR  TEMP.  AND  NINF' 

READ  *, TEMP, NINF 

PRINT  *, 'INPUT  VALUES  FOR  GAMMA  AND  b' 

READ  * , G , B 

IM=1  !  Number  of  equations 

Y(l)  =  1,00  !  Initial  condition  for  wl  at  r''=XP. 

Y(2)  =  0  I  Initial  condition  for  w2  at  r~=XP  (if  nec) 

PRINT  * 

PRINT  *, 'INTERVAL  OF  r~  FOR  PRINTING  ?' 

READ  *,PI 


PRINT  *, 'INPUT  THE  STEP  SIZE  (delta-r")' 
READ  *,H 


NS  =  NINT(PI/H) 


PRINT  *, 'MAXIMUM  r~  TO  STOP  CALCULATION  ?' 
READ  *,XL 


PRINT  * , '  dr"  =  ' ,H 

P1  =  0 
XP=0 
HH=H/2 . 

KPLUS=1 . 

print  * , '  NS=  ' ,  NS 
PRINT  * 


LI  =  0  !  Line  no,  initialization 

PRINT*, 'LINE  r~  KPLUS' 

WRITE  (*,98)  LI,XP,  KPLUS 
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28  LI=LI+1 

DO  N=1,NS 

XB=XP  !  OlJ  time 

XP=XP+H  I  New  time 

XM=XB+HH  I  Midpoint  time 

C  Runge-Kutta  Scheme - 

J=1  !  This  is  part  1. 

DO  1  =  1, IM 

yA( I )=Y( I ) 

END  no 
XA=XB 

CALL  FUNCT(rK, J, YA,K,XA,EINF,ERO,Q,ETA,G,B) 

J=2  i  This  is  part  2. 

DO  I=1,IM 

YA( I )=Y( I )+EK( 1 , I )/2 
END  DO 
XA=XM 

CALL  FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B ) 

J=3  !  This  is  part  3. 

DO  I  =  1  ,  I M 

YA( I )=Y( I )+EK( 2, I )/2 
END  DO 
XA=XM 

CALL  FUNCT( EK,J,YA,H, XA, EINF, ERO, Q, ETA, G,B) 

J=4  !  This  is  part  4. 

DO  1=1, IM 

YA( I )=Y( I )+EK( 3 , I ) 

END  DO 
XA=XP 

CALL  FUNCT ( EK , J , YA , H , XA , E INF , ERO , Q , ETA , G , B ) 

DO  1=1, IM  !  4-th  order  Runge-Kutta  scheme 

Y(I)=Y(I)  +  (EK(1,I)+EK(2,I  )  *2  +  EK(  3 , 1  )  >•  2  +  EK  (  4  ,  I  )  ) /6 
END  DO 

C  - 


W=Y( 1  ) 

c  You  now  have  W(  r~  )  —  to  get  K+(  r~  )  you  must 

c  multiply  by  {  E(  r~  )  /  E(0)  ) 

E=( EINF/ERO ) *DEXP(  ( Q**2 .  )/(  ( XP+1 .  )  **2  .  )  ) 

KPLUS=W*E 

c  To  get  the  derivative  of  W  (  W' (  r~  )  ] ,  just  plug  the 

c  computed  values  of  W  back  in  the  original  Ist-order  ODE 


J=5 

yA( 1 )=w 
XA=XP 

CALL  FUNCT (EK, J,YA,H,XA, EINF, ERO, Q, ETA , G , B ) 
WPRM=EK( 5, 1 )/H 
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n  h  n  '  'n  n  n  n  n 


NOW,  to  get  d/dr"  [ K+ ]  from  d/dr"  [  K+/(E/ErO)  ]  ,  must 
perform  the  following  operation: 

KPLSPRM=(WPRM*{E**2. ) +KPLUS*E* ( -2 . * ( Q**2 . )/( (XP+1. )**3. ) ) )/E 


AND  NOW,  to  generate  the  nondimensionalized  n-curves 
ne(  r""  )  and  ni  (  r~  ) 
compute  the  following  equations 

NI1=( KPLUS/E) *( EINF/ERO ) 

CC2=( K*TEMP*EPSO )/( (CHRG**2. ) *NINF* ( B* *2 . ) ) 
S2A=(6.*(Q**2. ) )/( (XP+1. )**4. ) 

S2B=(4.*(Q**4.  )  )/(  (XP+1.  )**6.  ) 

S2C=(1./(XP+(R0/B) )  )*( (-2.*(Q**2. ) )/( (XP+1.  )**3. )  ) 
S2D=(l./( (XP+(R0/B) )**2. ) ) 

S2=S2A+S2B+S2C-S2D 

NI2=CC2*S2 

CC3=( EPS0*ER0 )/( CHRG*NINF*B ) 

S3=( (-2.*(Q**2. ) )/( (XP+1. )**3. ) )+(l./(XP+(R0/B) ) ) 
NI3=CC3*E*S3 

NI=. 5* (NI1+NI2+NI3 ) 

CC1=( EPS0*ER0 )/( CHRG*NINF*B ) 

NE1=CC1*E*( (-2.*(Q**2. ) )/( (XP+1. )**3. ) + ( 1 . / ( XP+ ( RO/B ) ) ) ) 
=NI-NE1 


To  keep  from  generating  unplottable  50,000  point  data  files, 
the  following  edits  out  some  of  the  data  points  depending 
on  'where'  they  occur 

XCHK=(XP/H) 

IF  (XCHK  .GT.  100000)  GOTO  72 
IF  (XCHK  .GT.  10000)  GOTO  73 
IF  (XCHK  .GT.  1000)  GOTO  74 
IF  (XCHK  .GT.  100)  GOTO  75 
GOTO  44 

72  P1=P1+1 

IF  (PI  .NE.  10000)  GOTO  88 
P1  =  0 
GOTO  44 

73  P1=P1+1 

IF  (PI  .NE.  1000)  GOTO  88 
Pl  =  0 
GOTO  44 

74  Pl=Pl+l 

IF  (PI  .NE.  100)  GOTO  88 
P1  =  0 
GOTO  44 

75  P1=P1+1 

IF  (PI  .NE.  10)  GOTO  88 
P1  =  0 

44  WRITE  (10,*)  XP,  NI 

WRITE  (11,*)  XP,  NE 
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88 


CONTINUE 


END  DO 

WRITE  (*,98)  LI,XP,  KPLUS 
98  FORMATdX,  I2,  FlO.6,  2X,  1P4E16.8) 

IF  (XP  .LT.  XL)  GOTO  28 

200  PRINT* 

PRINT*, 'TYPE  1  TO  CONTINUE,  OR  0  TO  STOP,' 
READ  *,K1 

IF(Kl.EQ.l)  GOTO  1 
PRINT* 

END 


SUBROUTINE  FUNCT ( EK , J , YA , H , XA , EINF , ERO , Q , ETA , G , B )  !  DEFINES  ODES 

DOUBLE  PRECISION  EK ( 0 : 4 , 0 : 1 0 ) , YA ( 0 : 10 ) , H , XA , PART2 { 0 : 4 , 0 : 1 0 ) 
DOUBLE  PRECISION  PART3 ( 0 : 4 , 0 : 1 0 ) , PART4 ( 0 : 4 , 0 : 1 0 ) 

REAL  R0,KK4,KK3 
R0=.0001 

KK4  =  ETA*G*(  ( EINF/ERO ) * *2 .  ) 

PART4 ( J , 1 )=KK4*( DEXP( (2.*(Q**2.))/( (XA+1 . ) **2 . ) ) )/(XA*B+R0 ) 
KK3=ETA*( (G/R0)+1. )*( ( EINF/ERO )* *2 . ) 

PART3( J,1)=KK3*(DEXP( (2.*Q**2, )/( (XA+1. )**2. ) )/( (XA+1. )**3. ) ) 
PART2( J,1)=(DEXP( (Q**2. )/( (XA+1. )**2. ) ) )*( ETA* EINF/ERO ) 

EK(  J , 1 )  =  ( 2 . *ETA-PART2 ( J, 1 ) *YA( 1 ) -PART3 ( J , 1 ) +PART4 ( J, 1 )  ) *H 
C  EK(J,2)=  the  second  ode,  if  nec. 

RETURN 

END 
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